#!/usr/bin/perl

use DBI;
use POSIX;
use Cwd;
use Switch;
use strict;

# get our command line arguments
my %reduction_parameters=&parse_command_line(@ARGV);
my @rpfits_files;
if ($reduction_parameters{'rpfits_files'}){
    @rpfits_files=@{$reduction_parameters{'rpfits_files'}};
}
my @reduction_dates;
if ($reduction_parameters{'reduction_dates'}){
    @reduction_dates=@{$reduction_parameters{'reduction_dates'}};
}
my @project_codes;
if ($reduction_parameters{'project_codes'}){
    @project_codes=@{$reduction_parameters{'project_codes'}};
}

# database schema
my %observation;
($observation{'source'},$observation{'ra'},$observation{'dec'},$observation{'n_if'},
 $observation{'if_freq'},$observation{'if_nchans'},$observation{'if_bw'},$observation{'time'},
 $observation{'arrayname'},$observation{'project_code'},$observation{'rpfits_file'},
 $observation{'cycle_time'},$observation{'offsource'},$observation{'onsource'},
 $observation{'scantype'},$observation{'calcode'},$observation{'stations'},
 $observation{'header_number'},$observation{'num_processed'},$observation{'processed_ids'},
 $observation{'id'},$observation{'az'},$observation{'el'},$observation{'attenmm'},
 $observation{'subrefl'},$observation{'corr_cfg'},$observation{'airtemp'},
 $observation{'chi'},$observation{'evector'},$observation{'jyperk'},
 $observation{'pntra'},$observation{'pntdec'},$observation{'relhumid'},
 $observation{'smonrms'},$observation{'wind'},$observation{'winddir'},
 $observation{'xtsys'},$observation{'ytsys'},$observation{'xyphase'},$observation{'xcaljy'},
 $observation{'ycaljy'},$observation{'xgtp'},$observation{'ygtp'},$observation{'xsdo'},
 $observation{'ysdo'},$observation{'pointcor'},$observation{'point_time'},
 $observation{'point_source'},$observation{'point_az'},$observation{'point_el'},
 $observation{'usable'})=
    (0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,
     31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50);

# database details
my $dsn="dbi:Pg:dbname=atca;host=mentok;port=5432";
my $user="postgres";
my $password="";

# make the logfile
if (!$reduction_parameters{'output_logfile'}){
    $reduction_parameters{'output_logfile'}="imaging_".
	$reduction_parameters{'imaging_source'}.".log";
}
my $logname_determined=0;
my $test_logfile=$reduction_parameters{'output_logfile'};
my $append_number=1;
while ($logname_determined==0){
    if (-e $test_logfile){
	$test_logfile=$reduction_parameters{'output_logfile'}.".".
	    $append_number;
	$append_number++;
    } else {
	$logname_determined=1;
    }
}
$reduction_parameters{'output_logfile'}=$test_logfile;
open(LOGFILE,">".$reduction_parameters{'output_logfile'});
# make the LOG filehandle hot to prevent buffering
my $ofh = select LOGFILE;
$| = 1;
select $ofh;

# Setup the database connection
my $dbh=DBI->connect($dsn,$user,$password,{ RaiseError => 0, AutoCommit => 1});

&message_output(0,"MESSAGE","imager starts");
# output the arguments we're using
my $callstring=$0;
for (my $i=0;$i<=$#ARGV;$i++){
    $callstring.=" ".$ARGV[$i];
}
&message_output(1,"MESSAGE",$callstring);

# check we don't have conflicting arguments
if (($#reduction_dates>=0)&&($#rpfits_files>=0)){
    &message_output(-1,"ERROR","date was specified along with a list of RPFITS files");
    exit;
}

my $many_source_files=0;
if (($reduction_parameters{'mosaic_image'}==1)&&
    ($reduction_parameters{'individual_mosaic'}==1)){
    $many_source_files=1;
}
if (!$reduction_parameters{'mosaic_image'}){
    $reduction_parameters{'individual_mosaic'}=0; # just to be safe
}

my $completed=0;
my $current_date=0;
my @completed_directories;

while ($completed==0){
    # are we getting data from a specific date?
    if ($#reduction_dates>=0){
	$reduction_parameters{'date'}=$reduction_dates[$current_date];
	undef @rpfits_files;
	&message_output(0,"MESSAGE","working on data from date ".
			$reduction_parameters{'date'});
    }
    if ($reduction_parameters{'date'}){
	# find the earliest occurrence of the imaging source on that date
	my $ers_query="SELECT time,if_freq FROM observation WHERE ";
	if ($reduction_parameters{'mosaic_image'}==1){
	    $ers_query.="contains_text(source,'".$reduction_parameters{'imaging_source'}.
		"')=true";
	} else {
	    $ers_query.="source='".$reduction_parameters{'imaging_source'}."'";
	}
	$ers_query.=" AND time>='".$reduction_parameters{'date'}." 00:00:00' AND time<='".
	    $reduction_parameters{'date'}." 23:59:59'";
	if ($reduction_parameters{'reduce_frequency'}){
	    $ers_query.=" AND ".$reduction_parameters{'reduce_frequency'}.
		"=ANY(if_freq)";
	}
	$ers_query.=" ORDER BY time";
	my $ers=&start_database_handle($ers_query);
	my $earliesttime;
	my @ifconfig;
	my $nresults=0;
	while (my @ers_response=$ers->fetchrow){
	    $earliesttime=$ers_response[0];
	    @ifconfig=&postgresql_array_get($ers_response[1]);
	    $nresults++;
	    last;
	}
	$ers->finish;
	if ($nresults==0){
	    # we have no data for this source/frequency on this date
	    &message_output(0,"ERROR","no data present for source ".
			    $reduction_parameters{'imaging_source'}." with ".
			    "frequency ".$reduction_parameters{'reduce_frequency'}.
			    " on date ".$reduction_parameters{'date'});
	    @reduction_dates=&remove_date($current_date,@reduction_dates);
	    goto loop;
	}
	# find the nearest events before the first instance of that source on that date
	my @bounding_times=&find_nearest_events($earliesttime,@ifconfig);
	
	# get a list of all the RPFITS files between those times
	my $rpd_query="SELECT rpfits_file FROM observation WHERE ";
	if ($bounding_times[0] ne ""){
	    $rpd_query.="time>='".$bounding_times[0]."'";
	} elsif ($bounding_times[1] ne ""){
	    $rpd_query.="time>=(TIMESTAMP '".$bounding_times[1]."' - ".
		"INTERVAL '1 day')";
	} else {
	    $rpd_query.="time>=(TIMESTAMP '".$reduction_parameters{'date'}.
		" 00:00:00' - INTERVAL '1 day')";
	}
	$rpd_query.=" AND ";
	if ($bounding_times[1] ne ""){
	    $rpd_query.="time<='".$bounding_times[1]."'";
	} elsif ($bounding_times[0] ne ""){
	    $rpd_query.="time<=(TIMESTAMP '".$bounding_times[0]."' + ".
		"INTERVAL '1 day')";
	} else {
	    $rpd_query.="time<=(TIMESTAMP '".$reduction_parameters{'date'}.
		" 23:59:59' + INTERVAL '1 day')";
	}
	if ($#project_codes>=0){
	    $rpd_query.=" AND (";
	    for (my $i=0;$i<=$#project_codes;$i++){
		if ($i!=0){
		    $rpd_query.=" OR ";
		}
		$rpd_query.="project_code='".$project_codes[$i]."'";
	    }
	    $rpd_query.=")";
	}
	$rpd_query.=" ORDER BY time";
	my $rpd=&start_database_handle($rpd_query);
	while (my @rpd_response=$rpd->fetchrow){
	    my $already_on=0;
	    for (my $i=0;$i<=$#rpfits_files;$i++){
		if ($rpfits_files[$i] eq $rpd_response[0]){
		    $already_on=1;
		    last;
		}
	    }
	    if ($already_on==0){
		push @rpfits_files,$rpd_response[0];
		&message_output(2,"LIST","RPFITS file = ".$rpd_response[0]);
	    }
	}
	$rpd->finish;
    }

    my $nsource=1;
    my @source_list;
    my @all_sources;
    my @all_ids;
    my @all_starttimes;
    my @all_endtimes;
    my @all_ifs;
    # which mode are we operating in?
    if ($#rpfits_files>=0){
	# user wants us to make an image from a supplied set of files
	# make a list of all the sources
	my $lst_query="SELECT source,ra,dec,time,time+totaltime(maxcycles(onsource),".
	    "maxcycles(offsource),cycle_time),if_freq,id FROM observation WHERE (";
	for (my $i=0;$i<=$#rpfits_files;$i++){
	    if ($i!=0){
		$lst_query.=" OR ";
	    }
	    $lst_query.="contains_text(rpfits_file,'".$rpfits_files[$i]."')=true";
	}
	$lst_query.=") ORDER BY time";
	my $lst=&start_database_handle($lst_query);
	while (my @lst_response=$lst->fetchrow){
	    if ($reduction_parameters{'make_list'}==1){
		&message_output(0,"LIST","source $nsource (".$lst_response[3]." -> ".
				$lst_response[4]."): ".$lst_response[0]." [".
				$lst_response[1]." ".$lst_response[2]."] @ ".
				$lst_response[5]."\n");
	    }
	    @source_list=&add_source($lst_response[0],@source_list);
	    push @all_sources,$lst_response[0];
	    push @all_ids,$lst_response[6];
	    push @all_starttimes,$lst_response[3];
	    push @all_endtimes,$lst_response[4];
	    push @all_ifs,$lst_response[5];
	    $nsource++;
	}
	$lst->finish;
	
	if ($reduction_parameters{'make_list'}==1){
	    exit;
	}
    }

    # do some consistency checks
    # have we been given a source to image?
    if (!$reduction_parameters{'imaging_source'}){
	&message_output(-1,"ERROR","no source was specified for imaging");
	exit;
    }
    # are the specified sources in the files?
    for (my $i=0;$i<=$#source_list;$i++){
	&message_output(2,"LIST","found source ".$source_list[$i]);
    }
    my $source_found=&source_present($reduction_parameters{'imaging_source'},
				     $reduction_parameters{'mosaic_image'},
				     @rpfits_files);
    if ($source_found==0){
	&message_output(-1,"ERROR","imaging source is not present in the RPFITS files");
	exit;
    }
    if ($reduction_parameters{'phase_calibrator'}){
	$source_found=&source_present($reduction_parameters{'phase_calibrator'},0,
				      @rpfits_files);
	if ($source_found==0){
	    &message_output(-1,"ERROR","phase calibrator source is not present in the ".
			    "RPFITS files");
	    exit;
	}
    }
    if ($reduction_parameters{'flux_calibrator'}){
	$source_found=&source_present($reduction_parameters{'flux_calibrator'},0,
				      @rpfits_files);
	if ($source_found==0){
	    &message_output(-1,"ERROR","flux calibrator source is not present in the ".
			    "RPFITS files");
	    exit;
	}
    }
    if ($reduction_parameters{'bandpass_calibrator'}){
	$source_found=&source_present($reduction_parameters{'bandpass_calibrator'},0,
				      @rpfits_files);
	if ($source_found==0){
	    &message_output(-1,"ERROR","bandpass calibrator source is not present in ".
			    "the RPFITS files");
	    exit;
	}
    }    
    
    # check that the specified phase, flux and bandpass calibrators are phase, flux
    # and bandpass calibrators
    if ($reduction_parameters{'phase_calibrator'}){
	my $pc_type=&calibrator_type($reduction_parameters{'phase_calibrator'});
	if ($pc_type eq ""){
	    &message_output(-1,"ERROR","phase calibrator is not a known calibrator");
	    exit;
	}
    }
    if ($reduction_parameters{'flux_calibrator'}){
	my $pc_type=&calibrator_type($reduction_parameters{'flux_calibrator'});
	if ($pc_type ne "f"){
	    &message_output(-1,"ERROR","flux calibrator is not a known flux calibrator");
	    exit;
	}
    }
    if ($reduction_parameters{'bandpass_calibrator'}){
	my $pc_type=&calibrator_type($reduction_parameters{'bandpass_calibrator'});
	if ($pc_type eq ""){
	    if (!$reduction_parameters{'force_bandpass'}){
		&message_output(-1,"ERROR","bandpass calibrator is not a known calibrator");
		exit;
	    } else {
		&message_output(0,"WARNING","bandpass calibrator ".
				$reduction_parameters{'bandpass_calibrator'}." is likely ".
				"unsuitable; proceeding due to user override");
	    }
	}
    }

    if ($reduction_parameters{'mosaic_image'}==1){
	my $mosaic_message="imaging mosaic field ";
	if ($reduction_parameters{'individual_mosaic'}==1){
	    $mosaic_message.="(using individual approach) ";
	}
	&message_output(0,"MESSAGE",$mosaic_message.
			$reduction_parameters{'imaging_source'});
    } else {
	&message_output(0,"MESSAGE","imaging source ".$reduction_parameters{'imaging_source'});
    }

    # what frequencies are we imaging?
    my @unique_imaging_ifs;
    for (my $i=0;$i<=$#all_sources;$i++){
	if (($all_sources[$i] eq $reduction_parameters{'imaging_source'})||
	    (($reduction_parameters{'mosaic_image'}==1)&&
	     ($all_sources[$i]=~/^$reduction_parameters{'imaging_source'}/))){
	    my $already_printed=0;
	    for (my $j=0;$j<=$#unique_imaging_ifs;$j++){
		if ($unique_imaging_ifs[$j] eq $all_ifs[$i]){
		    $already_printed=1;
		    last;
		}
	    }
	    if ($already_printed==0){
		push @unique_imaging_ifs,$all_ifs[$i];
		&message_output(1,"MESSAGE","imaging source found with frequency ".
				"configuration ".$all_ifs[$i]);
	    }
	}
    }
    if (!$reduction_parameters{'reduce_frequency'}){
	if ($#unique_imaging_ifs>0){
	    # too many IF configurations, what would we do?!
	    &message_output(-1,"ERROR","too many IF configurations; you must specify which ".
			    "to image");
	    exit;
	}
	# choose the frequency of the first IF as the reduction frequency
	my @if_freqs=&postgresql_array_get($unique_imaging_ifs[0]);
	$reduction_parameters{'reduce_frequency'}=$if_freqs[0];
	$reduction_parameters{'reduce_if'}=1;
    } else {
	# check that the frequency supplied is one of the available frequencies
	my $found_frequency=0;
	for (my $i=0;$i<=$#unique_imaging_ifs;$i++){
	    my @if_freqs=&postgresql_array_get($unique_imaging_ifs[$i]);
	    for (my $j=0;$j<=$#if_freqs;$j++){
		if ($reduction_parameters{'reduce_frequency'}==$if_freqs[$j]){
		    $reduction_parameters{'reduce_if'}=$j+1;
		    $reduction_parameters{'reduce_ifconfig'}=$unique_imaging_ifs[$i];
		    $found_frequency=1;
		    last;
		}
	    }
	    if ($found_frequency==1){
		last;
	    }
	}
	if ($found_frequency==0){
	    # didn't find that frequency
	    &message_output(-1,"ERROR","could not identify any observations with frequency ".
			    $reduction_parameters{'reduce_frequency'});
	    exit;
	}
    }
    &message_output(0,"MESSAGE","will image frequency ".$reduction_parameters{'reduce_frequency'}.
		    " MHz (IF #".$reduction_parameters{'reduce_if'}.")");
    # is the bandpass calibrator appropriate for this frequency
    if ($reduction_parameters{'bandpass_calibrator'}){
	my $isbandpass=&is_bandpass_calibrator($reduction_parameters{'bandpass_calibrator'},
					       $reduction_parameters{'reduce_frequency'});
	if ($isbandpass==0){
	    if (!$reduction_parameters{'force_bandpass'}){
		# no good
		&message_output(-1,"ERROR","cannot use ".
				$reduction_parameters{'bandpass_calibrator'}." as a bandpass ".
				"calibrator at ".$reduction_parameters{'reduce_frequency'}.
				" MHz");
		exit;
	    } else {
		&message_output(0,"WARNING","bandpass calibrator ".
				$reduction_parameters{'bandpass_calibrator'}." is likely ".
				"unsuitable; proceeding due to user override");
	    }
	} else {
	    &message_output(0,"MESSAGE","using ".$reduction_parameters{'bandpass_calibrator'}.
			    " as the bandpass calibrator");
	}
    }

    # we now know exactly what configurations we're using, so we go through
    # and remake our RPFITS file list, and the source list
    # find out the bounding event times for the imaging source
    &message_output(1,"MESSAGE","redetermining RPFITS files");
    my $first_source=-1;
    for (my $i=0;$i<=$#all_sources;$i++){
	if (($all_sources[$i] eq $reduction_parameters{'imaging_source'})||
	    (($reduction_parameters{'mosaic_image'}==1)&&
	     ($all_sources[$i]=~/^$reduction_parameters{'imaging_source'}/))){
	    $first_source=$i;
	    last;
	}
    }
    my @source_ifs=&postgresql_array_get($reduction_parameters{'reduce_ifconfig'});
    my @near_events=&find_nearest_events($all_starttimes[$first_source],
					 @source_ifs);
    undef @rpfits_files;
    my $rpd_query="SELECT rpfits_file FROM observation WHERE ";
    if ($near_events[0] ne ""){
	$rpd_query.="time>='".$near_events[0]."'";
    }
    if ($near_events[1] ne ""){
	if ($near_events[0] ne ""){
	    $rpd_query.=" AND ";
	}
	$rpd_query.="time<='".$near_events[1]."'";
    }
    if ($#project_codes>=0){
	$rpd_query.=" AND (";
	for (my $i=0;$i<=$#project_codes;$i++){
	    if ($i!=0){
		$rpd_query.=" OR ";
	    }
	    $rpd_query.="project_code='".$project_codes[$i]."'";
	}
	$rpd_query.=")";
    }
    $rpd_query.=" AND if_freq='".$reduction_parameters{'reduce_ifconfig'}."'";
    $rpd_query.=" ORDER BY time";
    my $rpd=&start_database_handle($rpd_query);
    while (my @rpd_response=$rpd->fetchrow){
	my $already_on=0;
	for (my $i=0;$i<=$#rpfits_files;$i++){
	    if ($rpfits_files[$i] eq $rpd_response[0]){
		$already_on=1;
		last;
	    }
	}
	if ($already_on==0){
	    push @rpfits_files,$rpd_response[0];
	    &message_output(2,"LIST","RPFITS file = ".$rpd_response[0]);
	}
    }
    $rpd->finish;

    # now update the source list
    undef @source_list;
    undef @all_sources;
    undef @all_ids;
    undef @all_starttimes;
    undef @all_endtimes;
    undef @all_ifs;
    my $lst_query="SELECT source,ra,dec,time,time+totaltime(maxcycles(onsource),".
	"maxcycles(offsource),cycle_time),if_freq,id FROM observation WHERE (";
    for (my $i=0;$i<=$#rpfits_files;$i++){
	if ($i!=0){
	    $lst_query.=" OR ";
	}
	$lst_query.="contains_text(rpfits_file,'".$rpfits_files[$i]."')=true";
    }
    $lst_query.=") AND if_freq='".$reduction_parameters{'reduce_ifconfig'}."'".
	" ORDER BY time";
    my $lst=&start_database_handle($lst_query);
    while (my @lst_response=$lst->fetchrow){
	if ($reduction_parameters{'make_list'}==1){
	    &message_output(0,"LIST","source $nsource (".$lst_response[3]." -> ".
			    $lst_response[4]."): ".$lst_response[0]." [".
			    $lst_response[1]." ".$lst_response[2]."] @ ".
			    $lst_response[5]."\n");
	}
	@source_list=&add_source($lst_response[0],@source_list);
	push @all_sources,$lst_response[0];
	push @all_ids,$lst_response[6];
	push @all_starttimes,$lst_response[3];
	push @all_endtimes,$lst_response[4];
	push @all_ifs,$lst_response[5];
	$nsource++;
    }
    $lst->finish;
    for (my $i=0;$i<=$#source_list;$i++){
	&message_output(2,"LIST",$source_list[$i]);
    }
    for (my $i=0;$i<=$#all_sources;$i++){
	if ((($all_sources[$i] eq $reduction_parameters{'imaging_source'})||
	     (($reduction_parameters{'mosaic_image'}==1)&&
	      ($all_sources[$i]=~/^$reduction_parameters{'imaging_source'}/)))&&
	    ($all_ifs[$i] eq $reduction_parameters{'reduce_ifconfig'})){
	    $first_source=$i;
	    last;
	}
    }

    # figure out the best sources to use if the user didn't specify them
    # the phase calibrator
    if (!$reduction_parameters{'phase_calibrator'}){
	&message_output(0,"MESSAGE","attempting to determine most suitable phase calibrator");
	my @before_calibrators;
	my @after_calibrators;
	my $in_mosaic=0;
	for (my $i=0;$i<=$#all_sources;$i++){
	    if ((($all_sources[$i] eq $reduction_parameters{'imaging_source'})||
		 (($reduction_parameters{'mosaic_image'}==1)&&
		  ($all_sources[$i]=~/^$reduction_parameters{'imaging_source'}/)))&&
		($all_ifs[$i] eq $reduction_parameters{'reduce_ifconfig'})){
		if ($in_mosaic==1){
		    # we only look at one of the mosaic fields each time to save
		    # a lot of time searching through the database
		    next;
		}
		if ($reduction_parameters{'mosaic_image'}==1){
		    $in_mosaic=1;
		}
		my ($cal_before,$cal_after)=&nearest_source($all_starttimes[$i],
							    $all_endtimes[$i],$all_sources[$i],
							    $reduction_parameters{'reduce_ifconfig'},
							    $reduction_parameters{'mosaic_image'});
		push @before_calibrators,$cal_before;
		push @after_calibrators,$cal_after;
	    } elsif ($reduction_parameters{'mosaic_image'}==1){
		$in_mosaic=0;
	    }
	}
	my @count_source;
	my @source_count;
	for (my $i=0;$i<=$#before_calibrators;$i++){
	    my $new_source=1;
	    for (my $j=0;$j<=$#count_source;$j++){
		if ($count_source[$j] eq $before_calibrators[$i]){
		    $new_source=0;
		    $source_count[$j]++;
		    last;
		}
	    }
	    if ($new_source==1){
		push @count_source,$before_calibrators[$i];
		push @source_count,1;
	    }
	}
	for (my $i=0;$i<=$#after_calibrators;$i++){
	    my $new_source=1;
	    for (my $j=0;$j<=$#count_source;$j++){
		if ($count_source[$j] eq $after_calibrators[$i]){
		    $new_source=0;
		    $source_count[$j]++;
		    last;
		}
	    }
	    if ($new_source==1){
		push @count_source,$after_calibrators[$i];
		push @source_count,1;
	    }
	}
	my $max_count=0;
	my $max_position=-1;
	
	for (my $i=0;$i<=$#count_source;$i++){
	    if ($source_count[$i]>$max_count){
		$max_count=$source_count[$i];
		$max_position=$i;
	    }
	}
	$reduction_parameters{'phase_calibrator'}=$count_source[$max_position];
	&message_output(0,"MESSAGE","recommend using ".$reduction_parameters{'phase_calibrator'}.
			" as the phase calibrator");
    } else {
	&message_output(0,"MESSAGE","using ".$reduction_parameters{'phase_calibrator'}.
			" as the phase calibrator");
    }

    # the flux calibrator
    if (!$reduction_parameters{'flux_calibrator'}){
	&message_output(0,"MESSAGE","attempting to determine most suitable flux calibrator");
	my @fluxcal_ids=&find_nearest_fluxcals($all_starttimes[$first_source],
					       @near_events,@source_ifs);
	if ($#fluxcal_ids>=0){
	    my @id_response=&get_id_info($fluxcal_ids[0]);
	    $reduction_parameters{'flux_calibrator'}=$id_response[$observation{'name'}];
	    &message_output(0,"MESSAGE","recommend using ".$reduction_parameters{'flux_calibrator'}.
			    " as the flux calibrator");
	} else {
	    if (!$reduction_parameters{'relax_flux_calibration'}){
		&message_output(-1,"ERROR","unable to find suitable flux calibrator");
		exit;
	    } else {
		&message_output(0,"WARNING","unable to find suitable flux calibrator");
	    }
	}
	# check that the flux calibrator is part of the observations
	$source_found=&source_present($reduction_parameters{'flux_calibrator'},
				      @rpfits_files);
	if (($source_found==0)&&(!$reduction_parameters{'relax_flux_calibration'})){
	    &message_output(0,"MESSAGE","flux calibrator cannot be obtained from supplied ".
			    "files");
	    for (my $i=0;$i<=$#fluxcal_ids;$i++){
		my @id_response=&get_id_info($fluxcal_ids[$i]);
		&message_output(0,"MESSAGE","require file ".$id_response[$observation{'rpfits_file'}]);
		push @rpfits_files,$id_response[$observation{'rpfits_file'}];
		# add the sources from that file
		my $lst_query="SELECT source,ra,dec,time,time+totaltime(maxcycles(onsource),".
		    "maxcycles(offsource),cycle_time),if_freq,id FROM observation WHERE ".
		    "contains_text(rpfits_file,'".$id_response[$observation{'rpfits_file'}]."')".
		    "=true ORDER BY time";
		my $lst=&start_database_handle($lst_query);
		while (my @lst_response=$lst->fetchrow){
		    @source_list=&add_source($lst_response[0],@source_list);
		    push @all_sources,$lst_response[0];
		    push @all_ids,$lst_response[6];
		    push @all_starttimes,$lst_response[3];
		    push @all_endtimes,$lst_response[4];
		    push @all_ifs,$lst_response[5];
		    $nsource++;
		}
		$lst->finish;

		$source_found=&source_present($reduction_parameters{'flux_calibrator'},
					      @rpfits_files);
		if ($source_found==1){
		    last;
		}
	    }
	}
    } else {
	if ($reduction_parameters{'flux_calibrator'} ne ""){
	    &message_output(0,"MESSAGE","using ".$reduction_parameters{'flux_calibrator'}.
			    " as the flux calibrator");
	} else {
	    if ($reduction_parameters{'relax_flux_calibration'}==1){
		&message_output(0,"MESSAGE","will not attempt flux calibration");
	    } else {
		&message_output(-1,"ERROR","cannot flux calibrate this observation");
		&reduction_dates=&remove_date($current_date,@reduction_dates);
		goto loop;
	    }
	}
    }

    # the bandpass calibrator
    if (!$reduction_parameters{'bandpass_calibrator'}){
	&message_output(0,"MESSAGE","attempting to determine most suitable bandpass calibrator");
	for (my $i=0;$i<=$#source_list;$i++){
	    my $isbandpass=&is_bandpass_calibrator($source_list[$i],
						   $reduction_parameters{'reduce_frequency'});
	    if ($isbandpass==1){
		# we've found a suitable bandpass calibrator
		$reduction_parameters{'bandpass_calibrator'}=$source_list[$i];
		last;
	    }
	}
	if ($reduction_parameters{'bandpass_calibrator'}){
	    &message_output(0,"MESSAGE","using ".$reduction_parameters{'bandpass_calibrator'}.
			    " as the bandpass calibrator");
	} else {
	    if ($reduction_parameters{'replace_bandpass_calibration'}==1){
		&message_output(0,"MESSAGE","will attempt bandpass calibration with ".
				"potentially unsuitable calibrator");
	    } else {
		&message_output(-1,"ERROR","unable to find a suitable bandpass calibrator");
		&reduction_dates=&remove_date($current_date,@reduction_dates);
		goto loop;
	    }
	}
    }
    # if we have no bandpass calibrator, and the user is OK with it, we make the
    # phase calibrator the bandpass calibrator
    if (($reduction_parameters{'bandpass_calibrator'} eq "")&&
	($reduction_parameters{'replace_bandpass_calibration'}==1)){
	&message_output(0,"MESSAGE","using phase calibrator ".$reduction_parameters{'phase_calibrator'}.
			" to determine bandpass solution");
	$reduction_parameters{'bandpass_calibrator'}=
	    $reduction_parameters{'phase_calibrator'};
    }
    # special case here: if we have no flux calibrator, and the user is OK
    # with that, we make the bandpass calibrator the "flux" calibrator as well
    if (($reduction_parameters{'flux_calibrator'} eq "")&&
	($reduction_parameters{'relax_flux_calibration'}==1)){
	$reduction_parameters{'flux_calibrator'}=$reduction_parameters{'bandpass_calibrator'};
    }

    # we now attempt to find the RPFITS files
    my @rpfits_required;
    for (my $i=0;$i<=$#rpfits_files;$i++){
	$rpfits_files[$i]=~s/^.*\/([^\/]*)$/$1/;
	my $thisdirectory=getcwd;
	if (!-e $thisdirectory."/".$rpfits_files[$i]){
	    push @rpfits_required,splice @rpfits_files,$i,1;
	    $i--;
	} else {
	    $rpfits_files[$i]=$thisdirectory."/".$rpfits_files[$i];
	}
    }
    # check which of the required files actually has the sources we
    # require
    for (my $i=0;$i<=$#rpfits_required;$i++){
	my $file_required=0;
	my $flr_query="SELECT source FROM observation WHERE contains_text(rpfits_file,'".
	    $rpfits_required[$i]."')=TRUE AND if_freq='".$reduction_parameters{'reduce_ifconfig'}."'";
	my $flr=&start_database_handle($flr_query);
	while (my @flr_response=$flr->fetchrow){
	    if (($flr_response[0] eq $reduction_parameters{'phase_calibrator'})||
		($flr_response[0] eq $reduction_parameters{'flux_calibrator'})||
		($flr_response[0] eq $reduction_parameters{'bandpass_calibrator'})||
		(($flr_response[0] eq $reduction_parameters{'imaging_source'})||
		 (($reduction_parameters{'mosaic_image'}==1)&&
		  ($flr_response[0]=~/^$reduction_parameters{'imaging_source'}/)))){
		# this file is required
		$file_required=1;
		last;
	    }
	}
	$flr->finish;
	if ($file_required==0){
	    &message_output(1,"MESSAGE","RPFITS file ".$rpfits_required[$i].
			    " is not required for reduction");
	    splice @rpfits_required,$i,1;
	    $i--;
	}
    }
    for (my $i=0;$i<=$#rpfits_required;$i++){
	&message_output(1,"MESSAGE: attempting to find RPFITS file ".
			$rpfits_required[$i]);
	my $thisdirectory=getcwd;
	$rpfits_required[$i]=&obtain_rpfits($rpfits_required[$i],$thisdirectory);
	if (-e $rpfits_required[$i]){
	    # the file was copied successfully
	    push @rpfits_files,splice @rpfits_required,$i,1;
	    $i--;
	}
    }
    # did we get all the files?
    if ($#rpfits_required>=0){
	&message_output(-1,"ERROR","unable to locate all required RPFITS files");
	exit;
    }

    # we have our datasets, we know which calibrators to use, so let's get to it
    # make a directory if we need to
    my $directory_changed=0;
    my $date_directory;
    if ($#reduction_dates>=0){
	$date_directory=$reduction_parameters{'date'};
	if (-d $date_directory){
	    # the directory we need to make already exists
	    if ($reduction_parameters{'no_overwrite_directory'}==1){
		# the user doesn't want us to overwrite this directory
		&message_output(1,"MESSAGE","directory ".$date_directory." already ".
				"exists: we will not overwrite this data");
		$current_date++;
		goto loop;
	    }
	} else {
	    my $makedir_command="mkdir ".$reduction_parameters{'date'};
	    &run_command($makedir_command);
	}
	chdir $date_directory;
	$directory_changed=1;
    }
    
    # MIRIAD will output all files in lowercase, so we convert all our source names
    # to lowercase
    %reduction_parameters=&lowercase_names(0,%reduction_parameters);

    # first, we load the data into MIRIAD
    &message_output(0,"MESSAGE","loading data into MIRIAD format");
    my $uvdata=&miriad_load_data($reduction_parameters{'reduce_frequency'},
				 $reduction_parameters{'reduce_if'},
				 $reduction_parameters{'no_data_loading'},@rpfits_files);
    
    # then we flag based on calibration events that we know about
    &message_output(0,"MESSAGE","flagging setup data");
    my ($earliest_source_time,$latest_source_time)=
	&source_timefind($reduction_parameters{'imaging_source'},
			 $reduction_parameters{'reduce_ifconfig'},@rpfits_files);
    &event_flag($uvdata,$reduction_parameters{'reduce_ifconfig'},
		$earliest_source_time,$latest_source_time,@rpfits_files);
    
    # then we split the data
    &message_output(0,"MESSAGE","splitting dataset");
    my $mosaic_flag=0;
    if (($reduction_parameters{'mosaic_image'}==1)&&
	($reduction_parameters{'individual_mosaic'}!=1)){
	$mosaic_flag=1;
    }
    &miriad_split_data($uvdata,$reduction_parameters{'no_remove_data'},
		       $mosaic_flag);
    # check we have all the datasets at this point!
    my ($ce,$ck);
    if ($reduction_parameters{'individual_mosaic'}!=1){
	($ce,$ck)=&check_existence($reduction_parameters{'imaging_source'},
				   $reduction_parameters{'reduce_frequency'});
	if ($ce==1){
	    &message_output(0,"ERROR","dataset ".$ck." does not exist");
	    @reduction_dates=&remove_date($current_date,@reduction_dates);
	    goto loop;
	}
    } # otherwise the source names could be quite weird, and I can't think of a
      # way to check that they exist
    ($ce,$ck)=&check_existence($reduction_parameters{'phase_calibrator'},
			       $reduction_parameters{'reduce_frequency'});
    if ($ce==1){
	&message_output(0,"ERROR","dataset ".$ck." does not exist");
	@reduction_dates=&remove_date($current_date,@reduction_dates);
	goto loop;
    }
    ($ce,$ck)=&check_existence($reduction_parameters{'bandpass_calibrator'},
			       $reduction_parameters{'reduce_frequency'});
    if ($ce==1){
	&message_output(0,"ERROR","dataset ".$ck." does not exist");
	if ($reduction_parameters{'replace_bandpass_calibration'}==1){
	    # the user has said we can use the phase calibrator to determine
	    # the bandpass solution if we need to, and we need to right now
	    $reduction_parameters{'bandpass_calibrator'}=
		$reduction_parameters{'phase_calibrator'};
	    &message_output(0,"WARNING","will instead attempt to use the phase ".
			    "calibrator to determine bandpass solution");
	} else {
	    # we have to fail!
	    @reduction_dates=&remove_date($current_date,@reduction_dates);
	    goto loop;
	}
    }
    ($ce,$ck)=&check_existence($reduction_parameters{'flux_calibrator'},
			       $reduction_parameters{'reduce_frequency'});
    if ($ce==1){
	&message_output(0,"ERROR","dataset ".$ck." does not exist");
	if ($reduction_parameters{'relax_flux_calibration'}==1){
	    # the user doesn't care whether the observations are flux
	    # calibrated, so we continue on
	    $reduction_parameters{'flux_calibrator'}="";
	    &message_output(0,"WARNING","continuing without flux calibration");
	} else {
	    # we have to fail!
	    @reduction_dates=&remove_date($current_date,@reduction_dates);
	    goto loop;
	}
    }
    
    # determine the reference antenna
    $reduction_parameters{'reference_antenna'}=
	&determine_reference_antenna(@rpfits_files);
    &message_output(1,"MESSAGE","using antenna ".$reduction_parameters{'reference_antenna'}.
		    " as the reference");

    # first step in reduction is to calibrate the bandpass calibrator
    &message_output(0,"MESSAGE","determining bandpass solution");
    &message_output(1,"MESSAGE","calibrating");
    &miriad_calibrate($reduction_parameters{'bandpass_calibrator'},
		      $reduction_parameters{'reduce_frequency'},
		      $reduction_parameters{'reference_antenna'},
		      1,0,0);
    
    # and now we automatically flag bad data from this calibrator
    my $bandpass_flagged=0;
    if (!$reduction_parameters{'no_autoflag_bandpass'}){
	&message_output(1,"MESSAGE","flagging");
	&autoflag($reduction_parameters{'bandpass_calibrator'},
		  $reduction_parameters{'reduce_frequency'});
	$bandpass_flagged=1;
	
	# redo the calibration now with the flagged dataset
	&message_output(1,"MESSAGE","calibrating");
	&miriad_calibrate($reduction_parameters{'bandpass_calibrator'},
			  $reduction_parameters{'reduce_frequency'},
			  $reduction_parameters{'reference_antenna'},
			  1,0,0);
    }
    
    # next we copy the bandpass solution to the flux and phase calibrators
    &message_output(0,"MESSAGE","copying bandpass solution to flux and phase calibrators");
    if ($reduction_parameters{'bandpass_calibrator'} ne
	$reduction_parameters{'flux_calibrator'}){
	&miriad_copy_calibration($reduction_parameters{'bandpass_calibrator'},
				 $reduction_parameters{'flux_calibrator'},
				 $reduction_parameters{'reduce_frequency'});
    }
    
    if ($reduction_parameters{'bandpass_calibrator'} ne
	$reduction_parameters{'phase_calibrator'}){
	&miriad_copy_calibration($reduction_parameters{'bandpass_calibrator'},
				 $reduction_parameters{'phase_calibrator'},
				 $reduction_parameters{'reduce_frequency'});
    }
    
    # and calibrate both the flux and phase calibrators
    &message_output(0,"MESSAGE","calibrating flux and phase calibrators");
    &miriad_calibrate($reduction_parameters{'flux_calibrator'},
		      $reduction_parameters{'reduce_frequency'},
		      $reduction_parameters{'reference_antenna'},
		      0,1,1);
    &miriad_calibrate($reduction_parameters{'phase_calibrator'},
		      $reduction_parameters{'reduce_frequency'},
		      $reduction_parameters{'reference_antenna'},
		      0,1,0);

    # now flag these calibrators
    my $flux_flagged=0;
    if ((!$reduction_parameters{'no_autoflag_flux'})&&
	(($reduction_parameters{'bandpass_calibrator'} ne
	  $reduction_parameters{'flux_calibrator'})||
	 ($bandpass_flagged==0))){
	&message_output(1,"MESSAGE","automatically flagging flux calibrator");
	&autoflag($reduction_parameters{'flux_calibrator'},
		  $reduction_parameters{'reduce_frequency'});
	$flux_flagged=1;
    }
    my $phase_flagged=0;
    if ((!$reduction_parameters{'no_autoflag_phase'})&&
	(($reduction_parameters{'bandpass_calibrator'} ne
	  $reduction_parameters{'phase_calibrator'})||
	 ($bandpass_flagged==0))&&
	(($reduction_parameters{'flux_calibrator'} ne
	  $reduction_parameters{'phase_calibrator'})||
	 ($flux_flagged==0))){
	&message_output(1,"MESSAGE","automatically flagging phase calibrator");
	&autoflag($reduction_parameters{'phase_calibrator'},
		  $reduction_parameters{'reduce_frequency'});
	$phase_flagged=1;
    }

    # and recalibrate
    if ($flux_flagged==1){
	&message_output(1,"MESSAGE","calibrating flux calibrator");
	&miriad_calibrate($reduction_parameters{'flux_calibrator'},
			  $reduction_parameters{'reduce_frequency'},
			  $reduction_parameters{'reference_antenna'},
			  0,1,1);
    }
    if ($phase_flagged==1){
	&message_output(1,"MESSAGE","calibrating phase calibrator");
	&miriad_calibrate($reduction_parameters{'phase_calibrator'},
			  $reduction_parameters{'reduce_frequency'},
			  $reduction_parameters{'reference_antenna'},
			  0,1,0);
    }

    # now we set the correct flux scaling
    &message_output(0,"MESSAGE","correcting flux scale and bandpasses for phase ".
		    "calibrator");
    &miriad_fluxscale($reduction_parameters{'flux_calibrator'},
		      $reduction_parameters{'phase_calibrator'},
		      $reduction_parameters{'reduce_frequency'});
    

    # copy the final calibration to the imaging source
    &message_output(0,"MESSAGE","copying gain and bandpass calibration to imaging ".
		    "source");
    if ($reduction_parameters{'phase_calibrator'} ne
	$reduction_parameters{'imaging_source'}){
	&miriad_copy_calibration($reduction_parameters{'phase_calibrator'},
				 $reduction_parameters{'imaging_source'},
				 $reduction_parameters{'reduce_frequency'},
				 $many_source_files);
    }
    
    # and apply it
    &message_output(0,"MESSAGE","applying imaging source calibration solution");
    &miriad_apply_calibration($reduction_parameters{'imaging_source'},
			      $reduction_parameters{'reduce_frequency'},
			      $many_source_files);

    # change out of the directory if necessary
    push @completed_directories,$date_directory;
    if ($directory_changed==1){
	chdir "..";
    }

    # check whether we're done
    $current_date++;
  loop:   
    if ($current_date>$#reduction_dates){
	$completed=1;
    }

    %reduction_parameters=&lowercase_names(1,%reduction_parameters);

}

# make a combined dataset
%reduction_parameters=&lowercase_names(0,%reduction_parameters);
if ($#completed_directories>=0){
    &message_output(0,"MESSAGE","combining calibrated datasets");
    &miriad_combine_datasets($reduction_parameters{'imaging_source'},
			     $reduction_parameters{'reduce_frequency'},
			     $many_source_files,
			     @completed_directories);
    if ($reduction_parameters{'image_separately'}==1){
	# we'll make an image for each date's subset
	push @completed_directories,".";
    } else {
	# only making the combined image
	undef @completed_directories;
	push @completed_directories,".";
    }
}
%reduction_parameters=&lowercase_names(1,%reduction_parameters);

# calibration complete, so it's imaging time
for (my $i=0;$i<=$#completed_directories;$i++){    
    
    %reduction_parameters=&lowercase_names(0,%reduction_parameters);
    # reset the self calibration loop index
    $reduction_parameters{'self_calibration_loops'}=
	$reduction_parameters{'specified_self_calibration_loops'};
    
    # change directory if necessary
    my $directory_changed=0;
    if ($completed_directories[$i] ne "."){
	chdir $completed_directories[$i];
	$directory_changed=1;
    }
    
    my ($ra_axis,$dec_axis)=(-1,-1);
    my $imaging_completed=0;
    my $image_image;
    my %image_parameters;
    while ($imaging_completed==0){
	&message_output(0,"MESSAGE","creating dirty map");
	%image_parameters=
	    &miriad_make_image($reduction_parameters{'imaging_source'},$ra_axis,$dec_axis,
			       $reduction_parameters{'reduce_frequency'},
			       $reduction_parameters{'robust_value'},
			       $reduction_parameters{'mosaic_image'},$many_source_files,
			       %image_parameters);
	
	# now deconvolve
	&message_output(0,"MESSAGE","deconvolving image");
	$image_parameters{'no_mfclean'}=$reduction_parameters{'no_mfclean'};
	$image_parameters{'no_mossdi'}=$reduction_parameters{'no_mossdi'};
	%image_parameters=&miriad_deconvolve_image(%image_parameters);
	
	# and restore the image
	&message_output(0,"MESSAGE","restoring image");
	%image_parameters=&miriad_restore_image(%image_parameters);
	
	# now decide whether the image has been made correctly
	my @xpixels=@{$image_parameters{'xpixels'}};
	my @ypixels=@{$image_parameters{'ypixels'}};
	if ((($xpixels[0]!=-1)&&($ypixels[0]!=-1))||
	    ($reduction_parameters{'no_mfclean'}==1)||
	    ($reduction_parameters{'mosaic_image'}==1)){
	    $imaging_completed=1;
	}
	
	# make a measurement of the image size
	my $max_flux;
	%image_parameters=&miriad_image_size(%image_parameters);
	my @imagetripled=@{$image_parameters{'image_tripled'}};
	@xpixels=@{$image_parameters{'xpixels'}};
	@ypixels=@{$image_parameters{'ypixels'}};
	my @image=@{$image_parameters{'image'}};
	for (my $i=0;$i<=$#imagetripled;$i++){
	    if ($imagetripled[$i]==1){
		$xpixels[$i]/=3;
		$ypixels[$i]/=3;
		if (!$reduction_parameters{'no_image_crop'}){
		    # cut out the primary beam area
		    &message_output(0,"MESSAGE","cropping image");
		    &miriad_crop_image($image[$i],$xpixels[$i],$ypixels[$i]);
		}
	    }
	}
	@{$image_parameters{'xpixels'}}=@xpixels;
	@{$image_parameters{'ypixels'}}=@ypixels;
	
	# are we self-calibrating?
	if ($reduction_parameters{'self_calibration_loops'}>0){
	    $imaging_completed=0;
	    &message_output(0,"MESSAGE","performing self-calibration");
	    $image_parameters{'reference_antenna'}=$reduction_parameters{'reference_antenna'};
	    &miriad_self_calibrate(%image_parameters);
	    $reduction_parameters{'self_calibration_loops'}--;
	}	
	
	if ((($reduction_parameters{'mosaic_image'}==1)&&
	     ($image_parameters{'multifile'}!=1))||
	    ($reduction_parameters{'no_mfclean'}==1)){
	    # don't want to make triple-sized images
	    $ra_axis=-1;
	    $dec_axis=-1;
	}
    
	# we're out, so we've finished making the image!
	if ($image_parameters{'combined_image'}){
	    &message_output(0,"MESSAGE","imaging completed, image is in ".
			    $image_parameters{'combined_image'});
	} else {
	    my @image=@{$image_parameters{'image'}};
	    &message_output(0,"MESSAGE","imaging completed, image is in ".
			    $image[0]);
	}	    
	
	# change out of the directory if necessary
	if ($directory_changed==1){
	    chdir "..";
	}
	
	%reduction_parameters=&lowercase_names(1,%reduction_parameters);
	
    }
}

$dbh->disconnect;

# close the logfile
close(LOGFILE);

exit;

sub message_output {
    my ($debuglevel,$lineprefix,$message)=@_;

    # check whether this message goes on-screen
    if ($reduction_parameters{'verbosity'}>$debuglevel){
	print &timestamp().$lineprefix.": ".$message."\n";
    }
    # output the message to the logfile
    print LOGFILE &timestamp().$lineprefix.": ".$message."\n";

}

sub timestamp {
    my ($sec,$min,$hour,$mday,$mon,$year,$wday,
	$yday,$isdst)=localtime(time);
    my $timenow=sprintf "%4d-%02d-%02d %02d:%02d:%02d  ",
    $year+1900,$mon+1,$mday,$hour,$min,$sec;
    
    return $timenow;
}

sub postgresql_array_get {
    my ($arraystring)=@_;
    # turns a postgresql array string into a perl array
    $arraystring=~s/[\{\}]//g;
    my @els=split(/\,/,$arraystring);
    return @els;
}

sub find_nearest_events {
    my ($time,@ifs)=@_;

    my $time_difference="event_time-'".$time."'";
    my ($nearest_before,$nearest_after);
    my $evf_query="SELECT * FROM events WHERE (".
	"(event_type='acal' AND if_freq='".
	&array_to_text(@ifs)."')) AND EXTRACT(epoch FROM (".$time_difference."))<0".
	" ORDER BY abs(EXTRACT(epoch FROM (".$time_difference."))) LIMIT 1";
    my $evf=&start_database_handle($evf_query);

    while(my @evf_response=$evf->fetchrow){
	$nearest_before=$evf_response[1];
    }
    $evf->finish;

    my $evf_query="SELECT * FROM events WHERE (".
	"(event_type='acal' AND if_freq='".
	&array_to_text(@ifs)."')) AND EXTRACT(epoch FROM (".$time_difference."))>0".
	" ORDER BY abs(EXTRACT(epoch FROM (".$time_difference."))) LIMIT 1";
    my $evf=&start_database_handle($evf_query);

    while(my @evf_response=$evf->fetchrow){
	$nearest_after=$evf_response[1];
    }
    $evf->finish;

    return ($nearest_before,$nearest_after);
}

sub miriad_load_data {
    my ($frequency,$if_number,$no_reload_data,@rpfits_files)=@_;

    my $uv_dataset="loaded_data.uv";
    # remove any current dataset from the directory
    if ((-e $uv_dataset)&&($no_reload_data!=1)){
	my $remove_command="rm -rf ".$uv_dataset;
	&run_command($remove_command);
    }

    # ensure that there are no duplicates in the RPFITS file list
    my @new_rpfits_list;
    for (my $i=0;$i<=$#rpfits_files;$i++){
	my $found=0;
	for (my $j=0;$j<=$#new_rpfits_list;$j++){
	    if ($rpfits_files[$i] eq $new_rpfits_list[$j]){
		$found=1;
		last;
	    }
	}
	if ($found==0){
	    push @new_rpfits_list,$rpfits_files[$i];
	}
    }
    # and make sure the files are in time order
    @rpfits_files=sort @new_rpfits_list;
	
    if (($no_reload_data!=1)||(!-e $uv_dataset)){
	my $atlod_command="atlod in=";
	for (my $i=0;$i<=$#rpfits_files;$i++){
	    if ($i>0){
		$atlod_command.=",";
	    }
	    $atlod_command.=$rpfits_files[$i];
	}
	$atlod_command.=" out=".$uv_dataset;
	$atlod_command.=" options=birdie,rfiflag,noauto,xycorr,opcorr ";
	$atlod_command.="ifsel=".$if_number;
	
	my @loaddata_output=&run_command($atlod_command);
    
	&print_output_lines("MIRIAD",@loaddata_output);
    }

    return $uv_dataset;
}

sub run_command {
    my ($command)=@_;

    &message_output(2,"COMMAND",$command);

    # check we're not doing something unsafe
    if ($command=~/^rm\s/){
	if ($command=~/\s\*\s/){
	    # oh dear, we're going to rm the whole directory's
	    # contents, we'd better quit!
	    &message_output(2,"WARNING","unsafe command detected, aborting");
	    return;
	}
    }

    my @command_output;
    open(IN,"-|")||exec $command." 2>&1";
    while(<IN>){
	chomp;
	push @command_output,$_;
    }
    close(IN);

    return @command_output;
}

sub find_nearest_fluxcals {
    my ($time,$nearevent_before,$nearevent_after,@ifs)=@_;

    my $unixtime_reference=&time_to_unixtime($time);
    my @fluxcal_ids;
    my @fluxcal_sources;
    my @fluxcal_maxonsources;
    
    # find an observation with the correct IF configuration nearest to our
    # observation time that also falls between the known nearby event times
    my $endtime="time+totaltime(maxcycles(onsource),maxcycles(offsource),cycle_time)";
    my $fnd_query="SELECT observation.id,observation.source,observation.onsource FROM ".
	"observation,calibrators WHERE observation.source=calibrators.name AND ".
	"calibrators.type='f'";
    if ($nearevent_before ne ""){
	$fnd_query.=" AND (".$endtime.")>'".$nearevent_before."'";
    }
    if ($nearevent_after ne ""){
	$fnd_query.=" AND (".$endtime.")<'".$nearevent_after."'";
    }
    $fnd_query.=" AND if_freq='".&array_to_text(@ifs)."'".
	" ORDER BY abs(EXTRACT(epoch FROM ((".$endtime.")-'".$time."')))";
    my $fnd=&start_database_handle($fnd_query);
    while (my @fnd_response=$fnd->fetchrow){
	push @fluxcal_ids,$fnd_response[0];
	push @fluxcal_sources,$fnd_response[1];
	my @onsource_values=&postgresql_array_get($fnd_response[2]);
	my ($min_ons,$max_ons)=&minmax(@onsource_values);
	push @fluxcal_maxonsources,$max_ons;
    }
    $fnd->finish;

    # figure out which flux calibrator has the most time onsource
    my @fluxcal_found_source;
    my @fluxcal_found_onsource;
    for (my $i=0;$i<=$#fluxcal_sources;$i++){
	my $already_found=-1;
	for (my $j=0;$j<=$#fluxcal_found_source;$j++){
	    if ($fluxcal_found_source[$j] eq $fluxcal_sources[$i]){
		$already_found=$j;
		last;
	    }
	}
	if ($already_found==-1){
	    push @fluxcal_found_source,$fluxcal_sources[$i];
	    push @fluxcal_found_onsource,$fluxcal_maxonsources[$i];
	} else {
	    $fluxcal_found_onsource[$already_found]+=$fluxcal_maxonsources[$i];
	}
    }
    my ($best_source,$best_onsource);
    $best_source=$fluxcal_found_source[0];
    $best_onsource=$fluxcal_found_onsource[0];
    for (my $i=1;$i<=$#fluxcal_found_source;$i++){
	if ($fluxcal_found_onsource[$i]>$best_onsource){
	    $best_onsource=$fluxcal_found_onsource[$i];
	    $best_source=$fluxcal_found_source[$i];
	}
    }

    # now remove all ids that don't correspond to the best source
    for (my $i=0;$i<=$#fluxcal_ids;$i++){
	if ($fluxcal_sources[$i] ne $best_source){
	    splice @fluxcal_ids,$i,1;
	    splice @fluxcal_sources,$i,1;
	    splice @fluxcal_maxonsources,$i,1;
	    $i--;
	}
    }

    return (@fluxcal_ids);
}

sub minmax {
    # get the minimum and maximum values in an array
    my @array=@_;

    my $minval=$array[0];
    my $maxval=$array[0];
    for (my $i=0;$i<=$#array;$i++){
	if ($array[$i]<$minval){
	    $minval=$array[$i];
	}
	if ($array[$i]>$maxval){
	    $maxval=$array[$i];
	}
    }
    return ($minval,$maxval);
}

sub remove_date {
    my ($removeindex,@date_list)=@_;

    splice @date_list,$removeindex,1;

    return @date_list;
}

sub nearest_source {
    my ($start_time,$end_time,$source,$ifconfig,$is_mosaic)=@_;
    
    my ($source_before,$source_after)=("","");

    my $nrp_query="SELECT source FROM observation,calibrators WHERE EXTRACT(epoch FROM (".
	"time-'".$start_time."'))<0 AND ";
    if ($is_mosaic==1){
	$nrp_query.="contains_text(source,'".$source."')=false";
    } else {
	$nrp_query.="source!='".$source."'";
    }
    $nrp_query.=" AND name=source AND if_freq='".$ifconfig."' ".
	"ORDER BY EXTRACT(epoch FROM (time-'".$start_time."')) DESC LIMIT 1";
    my $nrp=&start_database_handle($nrp_query);
    while (my @nrp_response=$nrp->fetchrow){
	$source_before=$nrp_response[0];
	last;
    }
    $nrp->finish;

    my $nrl_query="SELECT source FROM observation,calibrators WHERE EXTRACT(epoch FROM (".
	"time-'".$end_time."'))>0 AND ";
    if ($is_mosaic==1){
	$nrl_query.="contains_text(source,'".$source."')=false";
    } else {
	$nrl_query.="source!='".$source."'";
    }
    $nrl_query.=" AND name=source AND if_freq='".$ifconfig."' ".
	"ORDER BY EXTRACT(epoch FROM (time-'".$end_time."')) LIMIT 1";
    my $nrl=&start_database_handle($nrl_query);
    while (my @nrl_response=$nrl->fetchrow){
	$source_after=$nrl_response[0];
	last;
    }
    $nrl->finish;

    return ($source_before,$source_after);
}

sub calibrator_type {
    my ($source)=@_;
    
    my $type="";
    my $chk_query="SELECT type FROM calibrators WHERE name='".$source."'";
    my $chk=&start_database_handle($chk_query);
    while (my @chk_response=$chk->fetchrow){
	$type=$chk_response[0];
    }
    $chk->finish;

    return $type;
}

sub is_bandpass_calibrator {
    my ($source,$frequency)=@_;

    my $wavelength=&which_wavelength($frequency);
    my $isbandpass=-1;
    my $chk_query="SELECT bandpass_".$wavelength." FROM calibrators WHERE name='".$source."'";
    my $chk=&start_database_handle($chk_query);
    while (my @chk_response=$chk->fetchrow){
	$isbandpass=$chk_response[0];
    }
    $chk->finish;

    return $isbandpass;
}

sub source_present {
    my ($source,$is_mosaic,@rpfits_files)=@_;

    my $is_present=0;
    if ($source eq ""){
	# just return now
	return $is_present;
    }
    
    my $chs_query="SELECT source FROM observation WHERE ";
    if ($is_mosaic==1){
	$chs_query.="contains_text(source,'".$source."')=true";
    } else {
	$chs_query.="source='".$source."'";
    }
    $chs_query.=" AND (";
    for (my $i=0;$i<=$#rpfits_files;$i++){
	$rpfits_files[$i]=~s/.*\/([^\/]*)$/$1/;
	if ($i!=0){
	    $chs_query.=" OR ";
	}
	$chs_query.="contains_text(rpfits_file,'".$rpfits_files[$i]."')=true";
    }
    $chs_query.=")";
    my $chs=&start_database_handle($chs_query);
    while (my @chs_response=$chs->fetchrow){
	if (($chs_response[0] eq $source)||
	    (($is_mosaic==1)&&
	     ($chs_response[0]=~/^$source/))){
	    $is_present=1;
	    last;
	}
    }
    $chs->finish;

    return $is_present;
}

sub source_timefind {
    my ($source,$ifspec,@rpfits_files)=@_;
    
    my $fnt_query="SELECT time FROM observation WHERE source='".$source.
	"' AND if_freq='".$ifspec."' AND (";
    for (my $i=0;$i<=$#rpfits_files;$i++){
	$rpfits_files[$i]=~s/.*\/([^\/]*)$/$1/;
	if ($i!=0){
	    $fnt_query.=" OR ";
	} 
	$fnt_query.="contains_text(rpfits_file,'".$rpfits_files[$i]."')=TRUE";
    }
    $fnt_query.=") ORDER BY time";
    my $fnt=&start_database_handle($fnt_query);
    my ($earliest_source_time,$latest_source_time);
    my $got_early=0;
    while (my @fnt_response=$fnt->fetchrow){
	if ($got_early==0){
	    $earliest_source_time=$fnt_response[0];
	    $got_early=1;
	}
	$latest_source_time=$fnt_response[0];
    }
    $fnt->finish;

    return ($earliest_source_time,$latest_source_time);
}

sub add_source {
    my ($new_source,@source_list)=@_;

    my $is_new=1;
    for (my $i=0;$i<=$#source_list;$i++){
	if ($source_list[$i] eq $new_source){
	    $is_new=0;
	    last;
	}
    }

    if ($is_new==1){
	push @source_list,$new_source;
    }

    return @source_list;

}

sub array_to_text {
    # turns an array into a postgresql compliant input
    my (@array)=@_;

    my $string="{";
    # check for letters, but we don't look for 'E' which may be used
    # in exponential numbers
    my $needsquotes=0;
    for (my $i=0;$i<=$#array;$i++){
	if ($array[$i]=~/[a-zA-DF-Z]/){
	    $needsquotes=1;
	}
    }
    for (my $i=0;$i<=$#array;$i++){
	if ($i>0){
	    $string.=",";
	}
	if ($needsquotes==1){
	    $string.="'";
	}
	if (!$array[$i]){
	    if ($needsquotes==1){
		$string.="NA";
	    } else {
		$string.="999";
	    }
	} else {
	    $string.=$array[$i];
	}
	if ($needsquotes==1){
	    $string.="'";
	}
    }
    $string.="}";
    return $string;
}

sub get_id_info {
    my ($get_id)=@_;

    my $get_query="SELECT * FROM observation WHERE id=".$get_id;
    my $get=&start_database_handle($get_query);
    my @get_response;
    while (@get_response=$get->fetchrow){
	last;
    }
    $get->finish;

    return @get_response;
}

sub time_to_unixtime {
    my ($ptime)=@_;
    
    $ptime=~/^(\d{4})\-(\d{2})\-(\d{2})\s(\d{2})\:(\d{2})\:(\d{2})/;
    my ($year,$month,$day,$hour,$minute,$second)=($1,$2,$3,$4,$5,$6);
    $year-=1900;
    $month-=1;

    my $utime=mktime($second,$minute,$hour,$day,$month,$year,0,0);

    return $utime;
}

sub which_wavelength {
    my ($frequency)=@_;

    my @highlow_21cm=(0,2000);
    my @highlow_13cm=(2000,4000);
    my @highlow_6cm=(4000,6500);
    my @highlow_3cm=(6500,10000);
    my @highlow_12mm=(15000,25000);
    my @highlow_7mm=(25000,55000);
    my @highlow_3mm=(80000,110000);

    my $wavelength;
    if (($frequency>=$highlow_21cm[0])&&
	($frequency<$highlow_21cm[1])){
	$wavelength="21cm";
    } elsif (($frequency>=$highlow_13cm[0])&&
	     ($frequency<$highlow_13cm[1])){
	$wavelength="13cm";
    } elsif (($frequency>=$highlow_6cm[0])&&
	     ($frequency<$highlow_6cm[1])){
	$wavelength="6cm";
    } elsif (($frequency>=$highlow_3cm[0])&&
	     ($frequency<$highlow_3cm[1])){
	$wavelength="3cm";
    } elsif (($frequency>=$highlow_12mm[0])&&
	     ($frequency<$highlow_12mm[1])){
	$wavelength="12mm";
    } elsif (($frequency>=$highlow_7mm[0])&&
	     ($frequency<$highlow_7mm[1])){
	$wavelength="7mm";
    } elsif (($frequency>=$highlow_3mm[0])&&
	     ($frequency<$highlow_3mm[1])){
	$wavelength="3mm";
    }

    return $wavelength;
}

sub miriad_split_data {
    my ($atlod_data,$keep_data,$is_mosaic)=@_;

    my $uvsplit_command="uvsplit vis=".$atlod_data." options=clobber";
    if ($is_mosaic==1){
	$uvsplit_command.=",mosaic";
    }
    my @uvsplit_output=&run_command($uvsplit_command);

    &print_output_lines("MIRIAD",@uvsplit_output);

    if ($keep_data!=1){
	# now remove the raw data, we won't be needing it anymore
	my $remove_command="rm -rf ".$atlod_data;
	&run_command($remove_command);
    }

}

sub print_output_lines {
    my ($line_prefix,@lines)=@_;
    
    for (my $i=0;$i<=$#lines;$i++){
	&message_output(2,$line_prefix,$lines[$i]);
    }
}

sub is_planet {
    my ($source)=@_;

    $source=~tr/A-Z/a-z/; # change to lowercase

    my $planet=0;
    my @planets=qw/mercury venus earth mars jupiter saturn uranus neptune/;

    for (my $i=0;$i<=$#planets;$i++){
	if ($source eq $planets[$i]){
	    $planet=1;
	    last;
	}
    }

    return $planet;
}

sub miriad_calibrate {
    my ($source,$frequency,$refant,$do_bandpass,$do_gains,$fluxcal)=@_;

    # if the source is a planet, we don't do any calibration
    if (&is_planet($source)==1){
	&message_output(1,"MESSAGE","calibration aborted; source is a planet");
	return;
    }

    if ($do_bandpass==1){
	# we solve for the bandpass using mfcal
	my $mfcal_command="mfcal vis=".$source.".".$frequency.
	    " refant=".$refant." interval=0.1";
	my @mfcal_output=&run_command($mfcal_command);
	&print_output_lines("MIRIAD",@mfcal_output);
    }
    if ($do_gains==1){
	# we solve for the gains using gpcal
	my $calibrated=0;
	my $use_qusolve=1;
	while ($calibrated==0){
	    $calibrated=1;
	    # delete the gains and leakage tables
	    my $delhd_command="delhd in=".$source.".".$frequency.
		"/gains";
	    my @delhd_output=&run_command($delhd_command);
	    &print_output_lines("MIRIAD",@delhd_output);
	    my $delhd_command="delhd in=".$source.".".$frequency.
		"/leakage";
	    my @delhd_output=&run_command($delhd_command);
	    &print_output_lines("MIRIAD",@delhd_output);

	    my $gpcal_command="gpcal vis=".$source.".".$frequency.
		" refant=".$refant." interval=0.1 options=xyvary";
	    if (($fluxcal==0)&&($use_qusolve==1)){
		$gpcal_command.=",qusolve";
	    }
	    my @gpcal_output=&run_command($gpcal_command);
	    &print_output_lines("MIRIAD",@gpcal_output);
	    # check that gpcal didn't fail
	    for (my $i=0;$i<=$#gpcal_output;$i++){
		if ($gpcal_output[$i]=~/^\#\#\#\sFatal\sError:.*degenerate$/){
		    $calibrated=0;
		    $use_qusolve=0;
		}
	    }
	}
    }
}

sub determine_reference_antenna {
    my (@rpfits_files)=@_;

    my $ref_query="SELECT onsource FROM observation WHERE (";
    for (my $i=0;$i<=$#rpfits_files;$i++){
	if ($i!=0){
	    $ref_query.=" OR ";
	}
	$ref_query.="contains_text(rpfits_file,'".$rpfits_files[$i]."')=true";
    }
    $ref_query.=")";
    my $ref=&start_database_handle($ref_query);
    my @total_onsource_cycles;
    my $initialised_count=0;
    while (my @ref_response=$ref->fetchrow){
	my @onsource_cycles=&postgresql_array_get($ref_response[0]);
	if ($initialised_count==0){
	    for (my $i=0;$i<=$#onsource_cycles;$i++){
		$total_onsource_cycles[$i]=$onsource_cycles[$i];
	    }
	    $initialised_count=1;
	} else {
	    for (my $i=0;$i<=$#onsource_cycles;$i++){
		$total_onsource_cycles[$i]+=$onsource_cycles[$i];
	    }
	}	    
    }
    $ref->finish;
    my ($minants,$maxants)=&minmax(@total_onsource_cycles);
    my @order_preference=(2,3,1,4,5); # never allow CA06 to be reference
    my $reference_antenna=-1;
    while($reference_antenna==-1){
	for (my $i=0;$i<=$#order_preference;$i++){
	    if ($total_onsource_cycles[$order_preference[$i]-1]==$maxants){
		# this antenna will work as the reference
		$reference_antenna=$order_preference[$i];
		last;
	    }
	}
	$maxants--;
    }
    
    return $reference_antenna;
}

sub get_indmos_file_list {
    my ($prefix,$suffix)=@_;

    my @file_list;
    # make a list of files in the current directory that match
    # the template prefix\d*.suffix
    # this is useful for determining the fields that have been
    # split from a mosaic while using the individual approach
    my $ls_command="ls -d ".$prefix."*.".$suffix;
    my @ls_output=&run_command($ls_command);
    for (my $i=0;$i<=$#ls_output;$i++){
	&message_output(4,"DEBUG",$ls_output[$i]);
	if (($ls_output[$i]=~/^$prefix[_]*\d+\.$suffix$/)&&
	    (-d $ls_output[$i])){
	    # this is what we're looking for
	    &message_output(3,"LIST","found mosaic field file ".
			    $ls_output[$i]);
	    push @file_list,$ls_output[$i];
	}
    }

    return @file_list;
}

sub miriad_copy_calibration {
    my ($calibrated_source,$uncalibrated_source,$frequency,$multifile)=@_;

    if ($multifile==0){
	my $gpcopy_command="gpcopy vis=".$calibrated_source.".".
	    $frequency." out=".$uncalibrated_source.".".$frequency;
	my @gpcopy_output=&run_command($gpcopy_command);
	&print_output_lines("MIRIAD",@gpcopy_output);
    } else {
	# get a list of the files
	my @file_list=&get_indmos_file_list($uncalibrated_source,$frequency);
	for (my $i=0;$i<=$#file_list;$i++){
	    my $gpcopy_command="gpcopy vis=".$calibrated_source.".".
		$frequency." out=".$file_list[$i];
	    my @gpcopy_output=&run_command($gpcopy_command);
	    &print_output_lines("MIRIAD",@gpcopy_output);
	}
    }
}

sub autoflag {
    my ($source,$frequency)=@_;

    my $max_flagging=30; # in percent, if the autoflagger flags more
                         # than this amount of data, it won't be allowed to
                         # apply the flagging
    my $absmax_flagging=95; # in percent, if the dataset is flagged this
                            # much after the autoflagger, it won't be applied

    my $miriad_mirflag_command="mirflag vis=".$source.".".$frequency.
	" int=600 options=rms,medtcal";
    &message_output(2,"MESSAGE","determining mirflag effect on dataset");
    my @mirflag_output=&run_command($miriad_mirflag_command.",noapply");
    &print_output_lines("MIRIAD",@mirflag_output);
    # check that the flagging wouldn't ruin the dataset
    my $apply_flagging=1;
    for (my $i=0;$i<=$#mirflag_output;$i++){
	if ($mirflag_output[$i]=~/Flagged\:\s+(.*?)\%\s+(.*)\%/){
	    my $flagged_before=$1;
	    my $flagged_after=$2;
	    if (($flagged_after>=$absmax_flagging)||
		(($flagged_after-$flagged_before)>=$max_flagging)){
		$apply_flagging=0;
		&message_output(2,"MESSAGE","mirflag operation not recommended");
	    }
	}
    }
    if ($apply_flagging==1){
	&message_output(2,"MESSAGE","applying mirflag flagging");
	@mirflag_output=&run_command($miriad_mirflag_command);
	&print_output_lines("MIRIAD",@mirflag_output);
    }

    $miriad_mirflag_command="mirflag vis=".$source.".".$frequency.
	" int=600";
    for (my $j=0;$j<2;$j++){
	my $mirflag_command=$miriad_mirflag_command;
	if ($j==0){
	    $mirflag_command.=" stokes=xx,xx,yy";
	} elsif ($j==1){
	    $mirflag_command.=" stokes=yy,xx,yy";
	}
	$mirflag_command.=" options=amp,medsed,short";
	@mirflag_output=&run_command($mirflag_command.",noapply");
	&print_output_lines("MIRIAD",@mirflag_output);
	# check that the flagging wouldn't ruin the dataset
	my $apply_flagging=1;
	for (my $i=0;$i<=$#mirflag_output;$i++){
	    if ($mirflag_output[$i]=~/Flagged\:\s+(.*?)\%\s+(.*)\%/){
		my $flagged_before=$1;
		my $flagged_after=$2;
		if (($flagged_after>=$absmax_flagging)||
		    (($flagged_after-$flagged_before)>=$max_flagging)){
		    $apply_flagging=0;
		    &message_output(2,"MESSAGE","mirflag operation not recommended");
		}
	    }
	}
	if ($apply_flagging==1){
	    &message_output(2,"MESSAGE","applying mirflag flagging");
	    @mirflag_output=&run_command($mirflag_command);
	    &print_output_lines("MIRIAD",@mirflag_output);
	}
    }
}

sub event_flag {
    my ($uvdata,$ifspec,$earlysource,$latesource,@rpfits_files)=@_;

    my @event_times;
    # remove any path information in the RPFITS files
    for (my $i=0;$i<=$#rpfits_files;$i++){
	$rpfits_files[$i]=~s/.*\/([^\/]*)$/$1/;
    }

    # find the earliest and latest times available in the RPFITS files
    my $tim_query="SELECT time FROM observation WHERE (";
    for (my $i=0;$i<=$#rpfits_files;$i++){
	if ($i!=0){
	    $tim_query.=" OR ";
	}
	$tim_query.="contains_text(rpfits_file,'".$rpfits_files[$i]."')=true";
    }
    $tim_query.=") AND if_freq='".$ifspec."' ORDER BY time LIMIT 1";
    my $tim=&start_database_handle($tim_query);
    my $earliest_time;
    while (my @tim_response=$tim->fetchrow){
	$earliest_time=$tim_response[0];
	last;
    }
    $tim->finish;
    &message_output(1,"MESSAGE","dataset starts at time ".$earliest_time);

    $tim_query="SELECT time+totaltime(maxcycles(onsource),maxcycles(offsource),".
	"cycle_time) FROM observation WHERE (";
    for (my $i=0;$i<=$#rpfits_files;$i++){
	if ($i!=0){
	    $tim_query.=" OR ";
	}
	$tim_query.="contains_text(rpfits_file,'".$rpfits_files[$i]."')=true";
    }
    $tim_query.=") AND if_freq='".$ifspec."' ORDER BY time DESC LIMIT 1";
    my $tim=&start_database_handle($tim_query);
    my $latest_time;
    while (my @tim_response=$tim->fetchrow){
	$latest_time=$tim_response[0];
	last;
    }
    $tim->finish;
    &message_output(1,"MESSAGE","dataset ends at time ".$latest_time);

    push @event_times,$earliest_time;
    # get the events that occurred during the files
    my $evs_query="SELECT * FROM events WHERE event_time>'".$earliest_time."' AND ".
	"event_time<'".$latest_time."' AND (event_type='dcal' OR ".
	"event_type='acal') AND if_freq='".$ifspec."' ORDER BY event_time";
    my $evs=&start_database_handle($evs_query);
    while(my @evs_response=$evs->fetchrow){
	my $event_time=$evs_response[1];
	push @event_times,$event_time;
    }
    $evs->finish;
    push @event_times,$latest_time;

    if ($#event_times==1){
	# there were no calibration events, so no flagging is required
	&message_output(1,"MESSAGE","no calibration events found, no flagging required");
	return;
    }

    # we now have a list of when each significant event occurred. we go through
    # this and determine the time difference between each event
    my $maxtimediff=0;
    my ($early_good_time,$late_good_time);
    for (my $i=0;$i<$#event_times;$i++){
	my $this_timediff=&time_difference($event_times[$i],
					   $event_times[$i+1]);
	if (&time_bigger($maxtimediff,$this_timediff)==1){
	    # ensure that the source time is greater than the earliest time
	    # and less than the latest time
	    if (((&time_earlier($event_times[$i],$earlysource)==1)||
		 (&time_earlier($event_times[$i],$latesource)==1))&&
		((&time_earlier($earlysource,$event_times[$i+1])==1)||
		 (&time_earlier($latesource,$event_times[$i+1])==1))){
		$maxtimediff=$this_timediff;
		$early_good_time=$event_times[$i];
		$late_good_time=$event_times[$i+1];
	    }
	}
    }
    # add a buffer to the earliest good time for safety
    my $cycle_time=&get_cycle_time(@rpfits_files);
    my $buffer_time=4*$cycle_time; # in seconds
    my $early_good_time=&time_add($early_good_time,$buffer_time);

    &message_output(1,"MESSAGE","biggest block of good data goes for ".$maxtimediff." between ".
		    $early_good_time." and ".$late_good_time);

    # now flag around this good time block
    if ($early_good_time ne $earliest_time){
	my $uvflag_command="uvflag vis=".$uvdata." flagval=flag \"select=time(".
	    &postgresql_to_miriad_time($earliest_time).",".
	    &postgresql_to_miriad_time($early_good_time).")\"";
	my @uvflag_response=&run_command($uvflag_command);
	&print_output_lines("MIRIAD",@uvflag_response);
    }

    if ($late_good_time ne $latest_time){
	my $uvflag_command="uvflag vis=".$uvdata." flagval=flag \"select=time(".
	    &postgresql_to_miriad_time($late_good_time).",".
	    &postgresql_to_miriad_time($latest_time).")\"";
	my @uvflag_response=&run_command($uvflag_command);
	&print_output_lines("MIRIAD",@uvflag_response);
    }	

}

sub time_difference {
    my ($early_time,$late_time)=@_;

    my $tdf_query="SELECT TIMESTAMP '".$late_time."' - TIMESTAMP '".$early_time."'";
    my $tdf=&start_database_handle($tdf_query);
    my $timediff;
    while (my @tdf_response=$tdf->fetchrow){
	$timediff=$tdf_response[0];
	last;
    }
    $tdf->finish;

    return $timediff;
}

sub time_earlier {
    my ($check_time,$reference_time)=@_;

    my $tme_query="SELECT TIMESTAMP '".$check_time."' < TIMESTAMP '".
	$reference_time."'";
    my $tme=&start_database_handle($tme_query);
    my $result;
    while (my @tme_response=$tme->fetchrow){
	$result=$tme_response[0];
	last;
    }
    $tme->finish;

    return $result;
}

sub time_bigger {
    my ($small_time,$big_time)=@_;

    my $tbg_query="SELECT INTERVAL '".$big_time."' > INTERVAL '".$small_time."'";
    my $tbg=&start_database_handle($tbg_query);
    my $result;
    while (my @tbg_response=$tbg->fetchrow){
	$result=$tbg_response[0];
	last;
    }
    $tbg->finish;

    return $result;
}

sub postgresql_to_miriad_time {
    my ($postgresql_time)=@_;
    
    $postgresql_time=~/^([^\-]*)\-([^\-]*)\-([^\-]*)\s([^\:]*)\:([^\:]*)\:(.*)$/;
    my ($year,$month,$day,$hour,$minute,$second)=($1,$2,$3,$4,$5,$6);
    
    my $miriad_time=sprintf "%4d%3s%02d:%02d:%02d:%02d",$year,
    &month_number_to_name($month),$day,$hour,$minute,$second;
    
    return $miriad_time;
}

sub month_number_to_name {
    my ($month_number)=@_;
    
    my @month_names=("JAN","FEB","MAR","APR","MAY","JUN",
		     "JUL","AUG","SEP","OCT","NOV","DEC");
    $month_number--;
    return $month_names[$month_number];
    
}

sub time_add {
    my ($reference_time,$interval)=@_;
    
    my $gtq_query="SELECT TIMESTAMP '".$reference_time."' + INTERVAL '".$interval." second'";
    my $gtq=&start_database_handle($gtq_query);

    my $final_time;
    while(my @times=$gtq->fetchrow){
	$final_time=$times[0];
    }
    $gtq->finish;

    return $final_time;
}

sub get_cycle_time {
    my (@rpfits_files)=@_;
    
    my $gct_query="SELECT cycle_time FROM observation WHERE rpfits_file='".
	$rpfits_files[0]."' LIMIT 1";
    my $gct=&start_database_handle($gct_query);
    my $cycletime;
    while (my @gct_response=$gct->fetchrow){
	$cycletime=$gct_response[0];
	last;
    }
    $gct->finish;

    return $cycletime;
}

sub miriad_self_calibrate {
#    my ($image_model,$image_map,$source,$frequency,$maximum_flux,$refant,$is_mosaic)=@_;
    my %image_parameters=@_;

    my @model=@{$image_parameters{'model'}};
    my @source_files=@{$image_parameters{'source_files'}};
    my @map=@{$image_parameters{'image_map'}};
    my @maximum_flux=@{$image_parameters{'max_flux'}};
    for (my $i=0;$i<=$#source_files;$i++){
	if (($image_parameters{'is_mosaic'}==1)&&
	    ($image_parameters{'multifile'}!=1)){
	    # self-calibrating a mosaic observation is a bit trickier than
	    # doing it for a single field!
	    
	    # start by splitting up the mosaic field into the individual
	    # chunks that will be self-calibrated
	    my $demos_template=$source_files[$i].".demos.";
	    if (-e $demos_template."1"){
		my $remove_command="rm -rf ".$demos_template."*";
		&run_command($remove_command);
	    }
	    my $demos_command="demos vis=".$source_files[$i].
		" map=".$map[$i]." out=".$demos_template.
		" options=detaper";
	    my @demos_output=&run_command($demos_command);
	    &print_output_lines("MIRIAD",@demos_output);
	    
	    # now self calibrate
	    my $selfcal_command="selfcal \"model=".$demos_template."*\" ".
		"vis=".$source_files[$i].
		" options=amplitude,mfs,mosaic interval=1 ".
		"clip=".(0.01*$maximum_flux[$i]);
	    my @selfcal_output=&run_command($selfcal_command);
	    &print_output_lines("MIRIAD",@selfcal_output);
	    
	    # and make sure the solution interval is large enough
	    my $puthd_command="puthd in=".$source_files[$i]."/interval".
		" value=1";
	    my @puthd_output=&run_command($puthd_command);
	    &print_output_lines("MIRIAD",@puthd_output);
	} else {
	    if ($image_parameters{'multifile'}!=1){
		my $gpscal_command="gpscal vis=".$source_files[$i].
		    " model=".$model[$i]." clip=".(0.01*$maximum_flux[$i]).
		    " interval=1 options=amplitude,mfs,xyvary ".
		    " refant=".$image_parameters{'reference_antenna'};
		my @gpscal_output=&run_command($gpscal_command);
		&print_output_lines("MIRIAD",@gpscal_output);
	    } else {
		# we still unfortunately have to use selfcal
		my $selfcal_command="selfcal model=".$model[$i]." ".
		    "vis=".$source_files[$i].
		    " options=amplitude,mfs interval=1 ".
		    "clip=".(0.01*$maximum_flux[$i]);
		my @selfcal_output=&run_command($selfcal_command);
		&print_output_lines("MIRIAD",@selfcal_output);
		# and make sure the solution interval is large enough
		my $puthd_command="puthd in=".$source_files[$i]."/interval".
		    " value=1";
		my @puthd_output=&run_command($puthd_command);
		&print_output_lines("MIRIAD",@puthd_output);
	    }
	}
    }

    # at the end of all this, if we're doing a mosaic field in individual
    # mode, then for fields with no calibratible emission, we won't have
    # done anything useful. the solution is to gpcopy the solutions from
    # each field and apply them to all the other fields
    for (my $i=0;$i<=$#source_files;$i++){
	for (my $j=0;$j<=$#source_files;$j++){
	    if ($i==$j){
		next;
	    }
	    my $gpcopy_command="gpcopy vis=".$source_files[$i].
		" out=".$source_files[$j]." mode=merge";
	    my @gpcopy_output=&run_command($gpcopy_command);
	    &print_output_lines("MIRIAD",@gpcopy_output);
	}
    }
}

sub miriad_fluxscale {
    my ($fluxcal,$phasecal,$frequency)=@_;

    # is the flux calibrator a planet?
    my $fluxplanet=&is_planet($fluxcal);

    # are the flux calibrator and the phase calibrator the same source?
    if ($fluxcal ne $phasecal){
	# set the gain scaling to be the same
	my $gpboot_command="gpboot vis=";
	if ($fluxplanet==0){
	    # flux calibrator should have the correct amplitude
	    $gpboot_command.=$phasecal.".".$frequency." cal=".$fluxcal.".".$frequency;
	} else {
	    # flux calibrator will not have been gain calibrated
	    $gpboot_command.=$fluxcal.".".$frequency." cal=".$phasecal.".".$frequency;
	}
	my @gpboot_output=&run_command($gpboot_command);
	&print_output_lines("MIRIAD",@gpboot_output);
    }

    # now correct the spectral slope
    my $mfboot_command="mfboot vis=".$fluxcal.".".$frequency;
    if ($fluxcal ne $phasecal){
	$mfboot_command.=",".$phasecal.".".$frequency;
    }
    $mfboot_command.=" \"select=source(".$fluxcal.")\"";
    my @mfboot_output=&run_command($mfboot_command);
    &print_output_lines("MIRIAD",@mfboot_output);
    
}

sub miriad_make_image {
    my ($imaging_source,$xpixels,$ypixels,$frequency,$robust_value,$is_mosaic,
	$multifile,%image_parameters)=@_;

#    my %image_parameters;
    $image_parameters{'is_mosaic'}=$is_mosaic;
    $image_parameters{'multifile'}=$multifile;
    $image_parameters{'source'}=$imaging_source;
    $image_parameters{'frequency'}=$frequency;
    undef @{$image_parameters{'image_tripled'}};
    if ($multifile==1){
	@{$image_parameters{'source_files'}}=
	    &get_indmos_file_list($imaging_source,$frequency);
    } else {
	push @{$image_parameters{'source_files'}},$imaging_source.".".$frequency;
    }

    my @source_files=@{$image_parameters{'source_files'}};
    if (($is_mosaic==1)&&($multifile==1)){
	# we need to find a reference position for all the images, and we take
	# this from the position of the first listed file
	my $uvlist_command="uvlist vis=".$source_files[0].
	    " options=source";
	my @uvlist_output=&run_command($uvlist_command);
	&print_output_lines("MIRIAD",@uvlist_output);
	for (my $i=0;$i<=$#uvlist_output;$i++){
	    if ($uvlist_output[$i]=~/^RA\:\s+(.*)$/){
		$image_parameters{'ref_ra'}=$1;
	    } elsif ($uvlist_output[$i]=~/^DEC\:\s+(.*)$/){
		$image_parameters{'ref_dec'}=$1;
	    }
	}
    }

    my @source_files=@{$image_parameters{'source_files'}};
    for (my $i=0;$i<=$#source_files;$i++){
	my $dataprefix=$source_files[$i];
	my $map=$dataprefix.".imap";
	my $beam=$dataprefix.".ibeam";
	# check if the datasets already exist, and delete them if they do
	if (-e $map){
	    my $remove_command="rm -rf ".$map;
	    &run_command($remove_command);
	}
	if (-e $beam){
	    my $remove_command="rm -rf ".$beam;
	    &run_command($remove_command);
	}

	my $invert_command="invert vis=".$dataprefix.
	    " map=".$map." beam=".$beam." stokes=i robust=".$robust_value;
	if (($xpixels!=-1)&&($ypixels!=-1)){
	    if ($image_parameters{'xpixels'}){
		my @pixels_x=@{$image_parameters{'xpixels'}};
		my @pixels_y=@{$image_parameters{'ypixels'}};
		$xpixels=$pixels_x[$i];
		$ypixels=$pixels_y[$i];
	    }
	    $xpixels*=3;
	    $ypixels*=3;
	    $invert_command.=" imsize=".$xpixels.",".$ypixels;
	    $invert_command.=" options=systemp,mfs,sdb";
	    push @{$image_parameters{'image_tripled'}},1;
	} else {
	    # this is probably the first time we've been asked to image, so
	    # we let MIRIAD do what it thinks is best
	    $invert_command.=" options=double,systemp,mfs";
	    if ($is_mosaic==1){
		$invert_command.=",mosaic";
		if ($multifile==1){
		    $invert_command.=" offset=".$image_parameters{'ref_ra'}.
			",".$image_parameters{'ref_dec'};
		}
	    }
	    push @{$image_parameters{'image_tripled'}},0;
	}

	my @invert_output=&run_command($invert_command);
	&print_output_lines("MIRIAD",@invert_output);
	
	# find the theoretical noise of the image
	my $theoretical_rms=-1;
	for (my $i=0;$i<=$#invert_output;$i++){
	    if ($invert_output[$i]=~/^Theoretical\srms\snoise\:\s(.*)$/){
		push @{$image_parameters{'theoretical_rms'}},$1;
	    }
	}

	# add the required parameters
	push @{$image_parameters{'image_map'}},$map;
	push @{$image_parameters{'image_beam'}},$beam;
	push @{$image_parameters{'image_prefix'}},$dataprefix;
	push @{$image_parameters{'xpixels'}},$xpixels;
	push @{$image_parameters{'ypixels'}},$ypixels;
    }

    # return with the name of the image, beam and theoretical noise, along
    # with the number of pixels in each image direction
#    return ($map,$beam,$xpixels,$ypixels,$imagetripled,@theoretical_rms);
    return %image_parameters;
}

sub miriad_apply_calibration {
    my ($source,$frequency,$multifile)=@_;

    if ($multifile==0){
	# use uvaver to apply the calibration tables to the visibility data
	my $uvaver_command="uvaver vis=".$source.".".$frequency.
	    " out=".$source.".".$frequency.".uvaver";
	my @uvaver_output=&run_command($uvaver_command);
	&print_output_lines("MIRIAD",@uvaver_output);
	
	# write the applied dataset over the unapplied one
	my $remove_command="rm -rf ".$source.".".$frequency;
	&run_command($remove_command);
	my $rename_command="mv ".$source.".".$frequency.".uvaver ".
	    $source.".".$frequency;
	&run_command($rename_command);
    } else {
	# get a list of the files
	my @file_list=&get_indmos_file_list($source,$frequency);
	for (my $i=0;$i<=$#file_list;$i++){
	    my $uvaver_command="uvaver vis=".$file_list[$i].
		" out=".$file_list[$i].".uvaver";
	    my @uvaver_output=&run_command($uvaver_command);
	    &print_output_lines("MIRIAD",@uvaver_output);
	
	    # write the applied dataset over the unapplied one
	    my $remove_command="rm -rf ".$file_list[$i];
	    &run_command($remove_command);
	    my $rename_command="mv ".$file_list[$i].".uvaver ".
		$file_list[$i];
	    &run_command($rename_command);
	}
    }

}

sub miriad_crop_image {
    my ($image,$xpixels,$ypixels)=@_;

    # cut out the central region xpixels x ypixels and overwrite the image
    my $new_image=$image.".new";
    my $halfx=floor($xpixels/2.0);
    my $halfy=floor($ypixels/2.0);
    my $imsub_command="imsub in=".$image." out=".$new_image.
	" \"region=relcenter,boxes(-".$halfx.",-".$halfy.",".$halfx.",".
	$halfy.")\"";
    my @imsub_output=&run_command($imsub_command);
    &print_output_lines("MIRIAD",@imsub_output);

    my $remove_command="rm -rf ".$image;
    &run_command($remove_command);

    my $move_command="mv ".$new_image." ".$image;
    &run_command($move_command);
}

sub miriad_deconvolve_image {
#    my ($source,$frequency,$map,$beam,$xpixels,$ypixels,$is_mosaic,$multifile,@rms)=@_;
    my %image_parameters=@_;
    my $max_attempts=5;
    my $niters=30;
    my $niters_fail_factor=2;

    my @source_files=@{$image_parameters{'source_files'}};
    my @dataprefix=@{$image_parameters{'image_prefix'}};
    my @xpixels=@{$image_parameters{'xpixels'}};
    my @ypixels=@{$image_parameters{'ypixels'}};
    my @map=@{$image_parameters{'image_map'}};
    my @beam=@{$image_parameters{'image_beam'}};
    my @rms=@{$image_parameters{'theoretical_rms'}};
    for (my $i=0;$i<=$#source_files;$i++){
#	my $dataprefix=$source.".".$frequency;
	my $model=$dataprefix[$i].".iclean";
	if (-e $model){
	    my $remove_command="rm -rf ".$model;
	    &run_command($remove_command);
	}
	my $clean_command;
	if ((($xpixels[$i]==-1)&&($ypixels[$i]==-1))||
	    ($image_parameters{'no_mfclean'}==1)){
	    my $converge_fail;
	    if (($image_parameters{'is_mosaic'}==1)&&
		($image_parameters{'multifile'}==0)){
		# we try to deconvolve a mosiac image
		for (my $nattempts=0;$nattempts<$max_attempts;$nattempts++){
		    if (-e $model){
			my $remove_command="rm -rf ".$model;
			&run_command($remove_command);
		    }
		    $clean_command="mosmem map=".$map[$i]." beam=".$beam[$i].
			" out=".$model." niters=".$niters;
		    my @clean_output=&run_command($clean_command);
		    &print_output_lines("MIRIAD",@clean_output);
		    # check for indication that deconvolution did not converge
		    $converge_fail=0;
		    for (my $i=0;$i<=$#clean_output;$i++){
			if ($clean_output[$i]=~/Failed to converge in NITERS iterations/){
			    $niters*=$niters_fail_factor;
			    $converge_fail=1;
			    last;
			}
		    }
		    if ($converge_fail==0){
			# we converged, stop now
			last;
		    }
		}
		if (($converge_fail==1)&&
		    ($image_parameters{'no_mossdi'}!=1)){
		    if (-e $model){
			my $remove_command="rm -rf ".$model;
			&run_command($remove_command);
		    }
		    # mosmem didn't work, so we try a mossdi
		    $clean_command="mossdi map=".$map[$i]." beam=".$beam[$i].
			" out=".$model." cutoff=".(5*$rms[$i])." niters=3000 options=positive";
		    my @clean_output=&run_command($clean_command);
		    &print_output_lines("MIRIAD",@clean_output);
		}
	    } else {
		# no idea about the image size, and invert was used to make a
		# normal dirty map, so the beam has only one plane, and we have
		# to use regular clean
		$clean_command="clean map=".$map[$i]." beam=".$beam[$i]." out=".$model.
		    " cutoff=".(5*$rms[$i])." speed=-0.5 niters=3000 options=positive";
		my @clean_output=&run_command($clean_command);
		&print_output_lines("MIRIAD",@clean_output);
	    }
	} else {
	    if ($image_parameters{'multifile'}==1){
		# we've made an individual mosaic!
		$clean_command="clean map=".$map[$i]." beam=".$beam[$i]." out=".$model.
		    " cutoff=".(5*$rms[$i])." speed=-0.5 niters=3000 options=positive";
		my @clean_output=&run_command($clean_command);
		&print_output_lines("MIRIAD",@clean_output);
	    } else {
		# invert was used to make a spectral dirty beam, so we use
		# mfclean
		$clean_command="mfclean map=".$map[$i]." beam=".$beam[$i]." out=".$model.
		    " cutoff=".(5*$rms[$i])." speed=-0.5 niters=3000 ";
		my $halfx=floor($xpixels[$i]/6.0);
		my $halfy=floor($ypixels[$i]/6.0);
		$clean_command.="\"region=relcenter,boxes(-".$halfx.",-".$halfy.",".
		    $halfx.",".$halfy.")\"";
	    }
	}
	if ($image_parameters{'is_mosaic'}!=1){
	    my @clean_output=&run_command($clean_command);
	    &print_output_lines("MIRIAD",@clean_output);
	}
	push @{$image_parameters{'model'}},$model;
    }

    # return the name of the model
#    return $model;
    return %image_parameters;
}

sub miriad_combine_datasets {
    my ($source,$frequency,$multifile,@directories)=@_;

    my @file_list;
    if ($multifile==1){
	# find the list of files from the first directory
	my $thisdir=cwd();
	chdir $directories[0];
	@file_list=&get_indmos_file_list($source,$frequency);
	chdir $thisdir;
    } else {
	push @file_list,$source.".".$frequency;
    }
    for (my $i=0;$i<=$#file_list;$i++){
	# remove any existing datasets
	my $output_dataset=$file_list[$i];
	if (-d $output_dataset){
	    my $remove_command="rm -rf $output_dataset";
	    &run_command($remove_command);
	}
	
	my $uvcat_command="uvcat vis=";
	for (my $j=0;$j<=$#directories;$j++){
	    if ($j!=0){
		$uvcat_command.=",";
	    }
	    $uvcat_command.=$directories[$j]."/".$file_list[$i];
	}
	$uvcat_command.=" out=".$output_dataset.
	    " options=unflagged";
	my @uvcat_output=&run_command($uvcat_command);
	&print_output_lines("MIRIAD",@uvcat_output);
    }
}

sub miriad_restore_image {
#    my ($source,$frequency,$map,$beam,$model)=@_;
    my %image_parameters=@_;

    my @source_files=@{$image_parameters{'source_files'}};
    my @dataprefix=@{$image_parameters{'image_prefix'}};
    my @xpixels=@{$image_parameters{'xpixels'}};
    my @ypixels=@{$image_parameters{'ypixels'}};
    my @map=@{$image_parameters{'image_map'}};
    my @beam=@{$image_parameters{'image_beam'}};
    my @model=@{$image_parameters{'model'}};
    my @rms=@{$image_parameters{'theoretical_rms'}};
    for (my $i=0;$i<=$#source_files;$i++){
	my $image=$dataprefix[$i].".irestor";
	if (-e $image){
	    my $remove_command="rm -rf ".$image;
	    &run_command($remove_command);
	}
	my $restor_command="restor model=".$model[$i]." beam=".$beam[$i].
	    " map=".$map[$i]." out=".$image;
	my @restor_output=&run_command($restor_command);
	&print_output_lines("MIRIAD",@restor_output);
	push @{$image_parameters{'image'}},$image;
    }

    # combine to form a big image if we're doing an individual-method mosaic
    if (($image_parameters{'is_mosaic'}==1)&&
	($image_parameters{'multifile'}==1)){
	my $final_image=$image_parameters{'source'}.".".
	    $image_parameters{'frequency'}.".combined";
	if (-e $final_image){
	    my $remove_command="rm -rf ".$final_image;
	    &run_command($remove_command);
	}
	my $linmos_command="linmos in=\"".$image_parameters{'source'}."*.irestor\"".
	    " out=".$final_image." options=taper";
	my @linmos_output=&run_command($linmos_command);
	&print_output_lines("MIRIAD",@linmos_output);
	$image_parameters{'combined_image'}=$final_image;
    }

    # return the size of the synthesised beam and the name of the image
#    return ($image,$ra_axis,$dec_axis);
    return %image_parameters;
}

sub miriad_image_size {
#    my ($image)=@_;
    my %image_parameters=@_;

#    my ($xpixels,$ypixels,$max_flux);
    my @xpixels;
    my @ypixels;
    my @max_flux;

    my @image=@{$image_parameters{'image'}};
    for (my $j=0;$j<=$#image;$j++){
	my $imlist_command="imlist in=".$image[$j]." options=statistics";
	my @imlist_output=&run_command($imlist_command);
	&print_output_lines("MIRIAD",@imlist_output);
	
	my $reached_planes=0;
	$max_flux[$j]=0;
	for (my $i=0;$i<=$#imlist_output;$i++){
	    my $line=$imlist_output[$i];
	    if ($line=~/Trc\=\((.*)\,(.*)\,(.*)\)/){
		$xpixels[$j]=$1;
		$ypixels[$j]=$2;
	    } elsif ($line=~/plane\s+Frequency/){
		$reached_planes=1;
	    } elsif ($reached_planes==1){
		$line=~s/^\s*//;
		my @els=split(/\s+/,$line);
		if ($els[3]>$max_flux[$j]){
		    $max_flux[$j]=$els[3];
		}
	    }
	}
    }

    @{$image_parameters{'xpixels'}}=@xpixels;
    @{$image_parameters{'ypixels'}}=@ypixels;
    @{$image_parameters{'max_flux'}}=@max_flux;

#    return ($xpixels,$ypixels,$max_flux);
    return %image_parameters;
}

sub primary_beamsize {
    my ($frequency)=@_;

    my $dishdiameter=22; # metres, the size of each ATCA antenna
    my $c=299792458; # m/s, the speed of light
    $frequency*=1e6; # Hz, the frequency of the observation

    my $wavelength=$c/$frequency; # m, the wavelength of the observation

    my $beamsize=$wavelength/$dishdiameter; # radians
}

sub obtain_rpfits {
    my ($rpfits_file,$work_directory)=@_;

    # try to find the RPFITS file
    my @search_paths=($work_directory,
		      "/DATA/ARCHIVE_1/ATCA/archive");
    my $found_file="";
    $rpfits_file=~s/^.*\/(.*)$/$1/;

    for (my $i=0;$i<=$#search_paths;$i++){
	$search_paths[$i]=~s/\/$//;
	my $search_file=$search_paths[$i]."/".$rpfits_file;
	if (-e $search_file){
	    $found_file=$search_file;
	    last;
	}
    }
    if ($found_file eq ""){
	# couldn't find this file!
	&message_output(1,"WARNING","unable to find RPFITS file ".$rpfits_file.
			" from local repository");

	# try to get this file from the ATOA
	$rpfits_file=~/(\d{4})\-(\d{2})\-(\d{2})\_(.*)/;
	my $rpfits_year=$1;
	my $rpfits_month=$2;
	my @rpfits_archive_search=("/data/ATOA_1/ATCA",
				   "/data/ATOA_2/ATCA");
	for (my $i=0;$i<=$#rpfits_archive_search;$i++){
	    my $rpfits_archive_path=$rpfits_archive_search[$i];
	    $rpfits_archive_path.="/".$rpfits_year."-".$rpfits_month;

	    my $copy_command="rsync -av ste616\@horus:".$rpfits_archive_path."/";
	    $copy_command.=$rpfits_file;
	    $copy_command.=" ".$work_directory."/";
	    &run_command($copy_command);
	    
	    $found_file=$work_directory."/".$rpfits_file;
	    if (!-e $found_file){
		$found_file="";
	    } else {
		last;
	    }
	}
    }

    if ($found_file eq ""){
	# still couldn't find the file
	&message_output(0,"WARNING","unable to find RPFITS file ".$rpfits_file.
			" from any repository");
	return "";
    }

    # return the location of the file we just extracted
    return $found_file;
    
}

sub print_help {
    my @arguments=@_;

    my $nshort=shift @arguments;
    my @short_arguments;
    for (my $i=0;$i<=$nshort;$i++){
	push @short_arguments,shift @arguments;
    }
    my @long_arguments;
    for (my $i=0;$i<=$nshort;$i++){
	push @long_arguments,shift @arguments;
    }
    my $noverrides=shift @arguments;
    my @short_overrides;
    for (my $i=0;$i<=$noverrides;$i++){
	push @short_overrides,shift @arguments;
    }
    my @long_overrides;
    for (my $i=0;$i<=$noverrides;$i++){
	push @long_overrides,shift @arguments;
    }
    my $nforces=shift @arguments;
    my @short_forces;
    for (my $i=0;$i<=$nforces;$i++){
	push @short_forces,shift @arguments;
    }
    my @long_forces;
    for (my $i=0;$i<=$nforces;$i++){
	push @long_forces,shift @arguments;
    }

    print "USAGE INFORMATION:\n";
    print " atca_db_automatic_imager_v4.pl [options] [RPFITS files]\n";
    print "\n where options are:\n";
    for (my $i=0;$i<=$nshort;$i++){
	print "  -".$short_arguments[$i]." / --".$long_arguments[$i]."\n";
	print "    ";
	switch ($i){
	    case 0 {
		print "increase the verbosity of the output";
	    }
	    case 1 {
		print "the name of the source to image";
	    }
	    case 2 {
		print "the name of the phase calibrator to use";
	    }
	    case 3 {
		print "the name of the flux calibrator to use";
	    }
	    case 4 {
		print "the name of the bandpass calibrator to use";
	    }
	    case 5 {
		print "print a list of sources in the RPFITS files";
	    }
	    case 6 {
		print "use data with this central frequency (in MHz) for the image";
	    }
	    case 7 {
		print "use data from this date";
	    }
	    case 8 {
		print "while imaging, set the 'robust' parameter to this value";
	    }
	    case 9 {
		print "use this number of self-calibration loops during imaging";
	    }
	    case 10 {
		print "create a logfile with this name";
	    }
	    case 11 {
		print "override one of the following default behaviours:\n";
		for (my $j=0;$j<=$noverrides;$j++){
		    print "     ".$short_overrides[$j]." / ".$long_overrides[$j]."\n";
		    print "      ";
		    switch ($j){
			case 0 {
			    print "don't redo the data loading steps";
			}
			case 1 {
			    print "don't delete data after using it";
			}
			case 2 {
			    print "don't attempt to autoflag the bandpass calibrator";
			}
			case 3 {
			    print "don't attempt to autoflag the phase calibrator";
			}
			case 4 {
			    print "don't attempt to autoflag the flux calibrator";
			}
			case 5 {
			    print "don't deconvolve using mfclean";
			}
			case 6 {
			    print "don't crop oversize images to the primary beam area";
			}
			case 7 {
			    print "don't re-reduce data from an already reduced date";
			}
			case 8 {
			    print "don't fail if the observations cannot be flux-calibrated";
			}
			case 9 {
			    print "use phase calibrator to determine bandpass solution if necessary";
			}
			case 10 {
			    print "prevent the use of mossdi even if mosmem fails to converge";
			}
		    }
		    print "\n";
		}
	    }
	    case 12 {
		print "use only data from this project";
	    }
	    case 13 {
		print "force the script to do one of the following things:\n";
		for (my $j=0;$j<=$nforces;$j++){
		    print "     ".$short_forces[$j]." / ".$long_forces[$j]."\n";
		    print "      ";
		    switch ($j){
			case 0 {
			    print "make an image for each date separately";
			}
			case 1 {
			    print "use the specified bandpass calibrator";
			}
			case 2 {
			    print "when mosaicing, use the 'individual' approach";
			}
		    }
		    print "\n";
		}
	    }
	    case 14 {
		print "display this usage information";
	    }
	    case 15 {
		print "the source to image is a mosaic observation";
	    }
	}
	print "\n";
    }    

}

sub parse_command_line {
    my @commandlineargs=@_;

    my %parameters;
    # set some defaults
    $parameters{'verbosity'}=0;
    $parameters{'robust_value'}=1;
    $parameters{'self_calibration_loops'}=0;
    if ($#commandlineargs>=0){
	$parameters{'print_help'}=0;
    } else {
	$parameters{'print_help'}=1;
    }

    my @short_options=("v","I","P","F","B","L","f","D","R","S","l","n","p","r","h","m");
    my @long_options=("verbose","image","phase","flux","bandpass","list",
		      "frequency","date","robust","selfcal","logfile",
		      "no","project","force","help","mosaic");
    my @overrides_short=("D","R","FB","FP","FF","MC","C","T","f","b","MS");
    my @overrides_long=("load","remove","autoflag-bandpass","autoflag-phase",
			"autoflag-flux","mfclean","crop","duplicatedate",
			"flux-calibration","bandpass-calibration","mossdi");
    my @force_short=("I","B","i");
    my @force_long=("imaging","bandpass","individual");
    for (my $i=0;$i<=$#commandlineargs;$i++){
	if ($commandlineargs[$i]=~/^\-(.*)/){
	    my $thisarg=$1;
	    for (my $j=0;$j<=$#short_options;$j++){
		if (($thisarg eq $short_options[$j])||
		    ($thisarg eq "-".$long_options[$j])){
		    # we recognise this option
		    switch ($j){
			case 0 {
			    # increase the verbosity
			    $parameters{'verbosity'}++;
			}
			case 1 {
			    # the name of the source to image
			    $parameters{'imaging_source'}=$commandlineargs[++$i];
			    $parameters{'specified_imaging_source'}=
				$parameters{'imaging_source'};
			}
			case 2 {
			    # the name of the phase calibrator
			    $parameters{'phase_calibrator'}=$commandlineargs[++$i];
			    $parameters{'specified_phase_calibrator'}=
				$parameters{'phase_calibrator'};
			}
			case 3 {
			    # the name of the flux calibrator
			    $parameters{'flux_calibrator'}=$commandlineargs[++$i];
			    $parameters{'specified_flux_calibrator'}=
				$parameters{'flux_calibrator'};
			}
			case 4 {
			    # the name of the bandpass calibrator
			    $parameters{'bandpass_calibrator'}=$commandlineargs[++$i];
			    $parameters{'specified_bandpass_calibrator'}=
				$parameters{'bandpass_calibrator'};
			}
			case 5 {
			    # the user wants a list of all the sources in the supplied
			    # files
			    $parameters{'make_list'}=1;
			}
			case 6 {
			    # reduce this frequency (in MHz)
			    $parameters{'reduce_frequency'}=$commandlineargs[++$i];
			}
			case 7 {
			    # look for files on this date
			    push @{$parameters{'reduction_dates'}},$commandlineargs[++$i];
			}
			case 8 {
			    # use this value for robust imaging
			    $parameters{'robust_value'}=$commandlineargs[++$i];
			}
			case 9 {
			    # use a self-calibration technique with this number of loops
			    $parameters{'self_calibration_loops'}=$commandlineargs[++$i];
			    $parameters{'specified_self_calibration_loops'}=
				$parameters{'self_calibration_loops'};
			}
			case 10 {
			    # output to a logfile
			    $parameters{'output_logfile'}=$commandlineargs[++$i];
			}
			case 11 {
			    # the user wants to override a default behaviour
			    my $this_override=$commandlineargs[++$i];
			    for (my $k=0;$k<=$#overrides_short;$k++){
				if (($this_override eq $overrides_short[$k])||
				    ($this_override eq $overrides_long[$k])){
				    switch ($k) {
					case 0 {
					    # don't redo loading steps
					    $parameters{'no_data_loading'}=1;
					}
					case 1 {
					    # don't remove files as we go along
					    $parameters{'no_remove_data'}=1;
					}
					case 2 {
					    # don't autoflag the bandpass data
					    $parameters{'no_autoflag_bandpass'}=1;
					}
					case 3 {
					    # don't autoflag the phase data
					    $parameters{'no_autoflag_phase'}=1;
					}
					case 4 {
					    # don't autoflag the flux data
					    $parameters{'no_autoflag_flux'}=1;
					}
					case 5 {
					    # don't use mfclean for deconvolution
					    $parameters{'no_mfclean'}=1;
					}
					case 6 {
					    # don't crop out the primary beam area in a
					    # triple-sized image
					    $parameters{'no_image_crop'}=1;
					}
					case 7 {
					    # if we find an existing date directory, we don't
					    # attempt to redo that date's reduction
					    $parameters{'no_overwrite_directory'}=1;
					}
					case 8 {
					    # don't worry about flux calibration
					    $parameters{'relax_flux_calibration'}=1;
					}
					case 9 {
					    # use the phase calibrator for the bandpass solution
					    # if there is no suitable bandpass calibrator
					    $parameters{'replace_bandpass_calibration'}=1;
					}
					case 10 {
					    # do not resort to using mossdi
					    $parameters{'no_mossdi'}=1;
					}
				    }
				}
			    }
			}
			case 12 {
			    # only use data from certain project codes
			    push @{$parameters{'project_codes'}},$commandlineargs[++$i];
			}
			case 13 {
			    # force a particular option
			    my $this_force=$commandlineargs[++$i];
			    for (my $k=0;$k<=$#force_short;$k++){
				if (($this_force eq $force_short[$k])||
				    ($this_force eq $force_long[$k])){
				    switch ($k) {
					case 0 {
					    # make an image for each subset date
					    $parameters{'image_separately'}=1;
					}
					case 1 {
					    # use whatever the user wants as the bandpass
					    # calibrator
					    $parameters{'force_bandpass'}=1;
					}
					case 2 {
					    # make an mosaic with the "individual" method
					    $parameters{'individual_mosaic'}=1;
					}
				    }
				}
			    }
			}
			case 14 {
			    # print the help information
			    $parameters{'print_help'}=1;
			}
			case 15 {
			    # the imaging source is actually a mosaic
			    $parameters{'mosaic_image'}=1;
			}
		    }
		}
	    }
	} else {
	    # must be an RPFITS file
	    push @{$parameters{'rpfits_files'}},$commandlineargs[$i];
	}
    }

    if ($parameters{'print_help'}==1){
	&print_help($#short_options,@short_options,@long_options,
		    $#overrides_short,@overrides_short,@overrides_long,
		    $#force_short,@force_short,@force_long);
	exit;
    }

    return %parameters;
}

sub lowercase_names {
    my ($mode,%parameters)=@_;

    if ($mode==0){
	# we're making the names lowercase and backing up the original names
	$parameters{'bandpass_calibrator'}=~tr/A-Z/a-z/;
	$parameters{'flux_calibrator'}=~tr/A-Z/a-z/;
	$parameters{'phase_calibrator'}=~tr/A-Z/a-z/;
	$parameters{'imaging_source'}=~tr/A-Z/a-z/;
    } elsif ($mode==1){
	# we're restoring the original names
	$parameters{'bandpass_calibrator'}=$parameters{'specified_bandpass_calibrator'};
	$parameters{'flux_calibrator'}=$parameters{'specified_flux_calibrator'};
	$parameters{'phase_calibrator'}=$parameters{'specified_phase_calibrator'};
	$parameters{'imaging_source'}=$parameters{'specified_imaging_source'};
    }

    return %parameters;
}    

sub start_database_handle {
    my ($query)=@_;

    my $works=0;

    my $query_handle;
    &message_output(2,"QUERY",$query);
    while($works==0){
	$query_handle=$dbh->prepare($query);
	my $handle_err=$query_handle->err;
	my $handle_errstring=$query_handle->errstr;
	$works=&manage_postgresql_errors($handle_err,$handle_errstring);
	if ($works==1){
	    $query_handle->execute;
	    my $handle_err=$query_handle->err;
	    my $handle_errstring=$query_handle->errstr;
	    $works=&manage_postgresql_errors($handle_err,$handle_errstring);
	}
    }

    return $query_handle;
}

sub manage_postgresql_errors {
    my ($errornum,$errorstring)=@_;

    &message_output(3,"PSQL",$errornum.": ".$errorstring);
    my $actioncode=1;

    if ($errornum eq ""){
	# no error
	return ($actioncode);
    }

    if (($errornum==7)&&
	($errorstring=~/server closed the connection unexpectedly/)){
	# we've lost database connection, so we try to reinit it
	&message_output(3,"PSQL","reconnecting to database");
	$dbh->disconnect;
	$dbh=DBI->connect($dsn,$user,$password,{ RaiseError => 0, AutoCommit => 1 });
	$actioncode=0;
    }

    # we return the new database handle
    return ($actioncode);
}

sub check_existence {
    my ($source,$frequency)=@_;

    my $returncode=0;
    if ($source eq ""){
	# no source
	return (0,"");
    }

    my $dataset=$source.".".$frequency;
    if (!-d $dataset){
	# this dataset doesn't exist
	$returncode=1;
    }

    return ($returncode,$dataset);
}
